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Abstract 


This report examines the feasibility of using a periodic 
dielectric layer, composed of alternating bars having dielectric 
constants and e 2 , am a frequency selective sub-reflector in 
order to permit feed separation in large aperture reflecting 
antenna systems. For oblique incidence, it is found that total 
transmission and total reflection can be obtained at different 
frequencies for proper choice of e 1# e 2 and the geometric 
parameters. The frequencies of total reflection and transmission 
can be estimated form wave phenomena occurring in a layer of 
uniform dielectric constant equal to the average for the periodic 
layers. About some of the frequencies of total transmission, the 
bandwidth for 90* transmission is found to be 40*. However, the 
bandwidth for 90* reflection is always found to be much narrower; 
the greatest value found being 2.5*. 



I. 


Introduction 


Separation of the feed structures for different frequency bands in 
large reflecting antennas has been achieved using sub-reflectors whose 
transmission and reflection coefficients are frequency dependent. For 
frequencies in one band, the sub-reflector acts as a perfect reflector, 
while for frequencies in another distinct band the sub-reflector is 
transparent to the radiation, thus permitting direct illumination of the 
main reflector by the feed. To date, periodic arrays of conducting 
plates, or apertures in a conductive screen, have been used as the fre- 
quency selective surface [1-4]. Typically, the conductors are placed on 
a dielectric layer that provides mechanical support. 

Use of a dielectric layer with periodically varying dielectric 
constant has been suggested as an alternative way to obtain a frequency 
selective surface. As considered here, the layer is composed of alter- 
nating strips of two materials having different dielectric constants, as 
shown in Figure 1. At mm frequencies, such dielectric layers offer the 
advantage of low absorption loss as compared to metallic screens. Since 
the layer thickness is on the order of a wavelength, the amount of mate- 
rial required would not be excessive at these high frequencies. 

This report describes a theoretical study of frequency selective 
reflection and transmission at dielectric layers of the type shown in 
Figure 1. Because this effort was intended as a limited feasibility 
study, we consider the case when the plane is incident perpendicular to 
the strips, and assume the electric field to be polarized along the 
strips, as in Figure 1. The layer is found to exhibit the desired fre- 
quency selective properties. It is also found that the approximate 
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Figure 1. Frequency selective surface composed of a periodic array 
of dielectric strips. 


locations of reflection and transmission bands can be predicted from wave 
properties of a layer having uniform dielectric constant equal to the 
average of that in the periodic layer. 

The significant wave properties are discussed qualitatively in 
Section II. In Section III, the necessary mathematical analysis is 
carried out to permit numerical evaluation of the reflection and trans- 
mission coefficients. Numerical results are presented in Section IV. 
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II . Wave Mechanisms for Frequency Selective Behavior 

In this section, we describe the wave phenomena that can provide 
frequency selective reflection and transmission at a periodic dielectric 
layer. The description given here is intended to clarify the nature of 
the subsequent analysis, and to provide a context for discussing the 
numerical results that have been obtained. 

Consider first a dielectric that is periodic along x but infinite 
along z, as shown in Figure 2. For two dimensional propagation in the 
(x,z) plane, the dielectric will support an infinite set of modes with 
different wavenumbers, « n along z, but each having the same Bloch wave- 
number kjj along x [5,6]. At low frequencies, only the lowest n = 0 mode 
will have a real wavenumber * 0 , while all other modes will be cut off (« n 
imaginary or complex). At somewhat higher frequencies, the n = -1 mode 
will also propagate (k.j real), while higher modes remain cut off. 

Further increase in frequency will result in more propagating modes. 

Consider now a semi-infinite, periodic dielectric illuminated by a 
plan wave incident from vacuum, as shown in Figure 3. The incident wave 
will excite all of the modes of the periodic structure. At a low enough 
frequency f ^ , only the n * 0 mode will propagate along z, as suggested in 
Figure 3a. Higher modes will decay exponentially away from the surface z 
= 0. However, at a higher frequency f 2 , the n = 0 and n = -1 modes can 
propagate. If and e 2 are not close to the dielectric constant of free 
space, two modes can propagate in the dielectric even for the periodicity 
d small enough compared to the free space wavelength Aq so that no grat- 
ing lobes are present in the field reflected into the region z < 0. In 
this case, the reflected field propagates only at the specular angle, as 
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Figure 2. Periodic dielectric of infinite extent along z showing two of 

the infinite set of Bloch waves. The dielectric and fields are 
assumed to have no variation along y (out of the plane of the 
paper) . 



lane wave for the cases of: a) low frequency f- at which only the n=0 
loch wave propagates along z; and b) higher frequency f_ at which the 
=0 and n=-l Bloch waves propagate along z. 



indicated in Figure 3b. If the periodic dielectric is of finite thick- 
ness h, as in Figure 4, the propagating modes will excite a plane wave in 
the vacuum region z > h below the dielectric. 

At low frequencies f lt only one mode propagates along z with real 
wavenumber k q , so that the layer acts approximately as if it had a uni- 
form dielectric constant equal to the average of that for the periodic 
layer. Thus the transmission properties will be similar to those of a 
uniform layer. In particular, the reflection coefficient will vanish at 
about the frequency for which * 0 h » n. In this case, the dielectric 
acts as a half-wave window, and there will be total transmission of the 
incident plane wave, as suggested in Figure 4a. 

At a higher frequency f 2 , both the n = 0 and n = -1 modes will 
propagate along z. These modes are excited at the top surface of the 
layer by the incident wave. When each mode reaches the bottom surface, 
it excites both modes traveling back to the top, as well as a transmitted 
plane wave in the air. Because of the phase-matching conditions at the 
top and bottom surfaces, the transmitted plane wave in the air propagates 
in the same direction as the incident plane wave. 

The modes in the layer that are traveling back towards the top 
surface excite a reflected plane wave in the air above the layer, as well 
as being scattered back into the layer, as suggested in Figure 4b. 
Repetition of this scattering process establishes the total field in the 
layer, and the total reflected and transmitted plane waves in the air. 

For some frequency f 2 , the phases of the two modes in the layer will be 
such as to add destructively in producing the transmitted plane wave, 
while constructively adding for the reflected plane wave, thereby produc- 
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periodically varying dielectric layer by: a) taking h=TT/ic , so that the 
layer acts as a half-wave window at the low frequency f- : and b) chosin 
dielectric and geometric parameters to achieve cancellation of the 
transmitted wave at a higher frequency f Q . 



ing the desired frequency selective property. However, at other frequen- 
cies in the range where the n * 0 and n ■ -1 modes propagate along z, the 
phases of these two modes may be such as to add for the transmitted plane 
wave and cancel for the reflected wave. Thus it is possible to have 
multiple frequencies for which total reflection and total transmission 
take place. 

The foregoing behavior is known in other related geometries to be 
associated with the excitation of waves guided along the layer [7,8]. To 
understand the connection with the guided waves, consider a layer of 
uniform dielectric constant equal to an average permittivity of the 
periodic layer defined by 

e a = < s l d l + e 2 d 2 ) / d • U) 

This layer will support guided waves whose fields vary sinusoidally in 
the layer, and decay away from the layer in the air [9]. For TE guided 
wave modes, the normalized wavenumber <3 g h with g=0,l,2 is plotted along 
the horizontal axis in Figure 5 versus the normalized free space wave- 
number k Q h = a>h/c , which is plotted along the vertical axis for e a = 2 . 

Because /3 g is greater than « 0 , there are no angles of incidence 9 at 
which a plane wave can satisfy the phase-match condition k Q sin0 = 0 g for 
direct excitation of the waveguide modes. However, excitation is pos- 
sible if the dielectric constant of the layer is a periodic function of 
x. In this case, the fields of each waveguide mode will consist of a 
series of space harmonics, one of which has fields that are very similar 
to those of a mode in a uniform layer having the same average dielectric 
constant « a . This space harmonic is designated m ■ 0, and has wavenumber 
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Figure 5. Dispersion curves for the g=0,l,2,3 modes guided by a uniform layer 
of relative dielectric constant 2. 




£g 0 along x that is close to /3g for the uniform layer. Other space har- 
monics have wavenumbers along x given by 

$gq - |3g 0 ♦ q2z/d, (q - ±1 , ±2 . (2) 

While 0gQ * (3g > k Q , it is possible to choose the periodicity d such 
that I0g f _il < k Q for m = -1. In this case, a plane wave incident at an 
angle 8 = sin -1 ( I /3g, I /k Q ) will couple to the space harmonic, and 
through it excite the waveguide mode. Once excited, this mode will re- 
radiate plane waves into the air regions above and below the layer 
through the same space harmonic. The process of excitation and re-radia- 
tion is depicted in Figure 6a for the case when k Q + > 2k / d > $g. 

This condition is sufficient to guarantee that only one space harmonic 
will give rise to a plane wave propagating away from the layer, and also 
implies that the plane wave propagates backward with respect to the 
direction of the waveguide mode. Guided waves that radiate some of their 
energy as they propagate are known as leaky waves [10]. 

The same physical processes hold if the incident wave is from the 
left, as shown in Figure 6b, except that the direction of propagation 
along x is reversed for the waveguide mode. The re-radiated plane wave 
above the layer adds to the reflected plane wave generated directly at 
the top surface of the layer to give the total reflected field. When the 
two components are in phase, strong reflections take place. However, 
when they are out of phase the reflected field is small and strong 
transmission occurs. Because the phases are frequency dependent, the 
overall reflection can have the desired frequency selective behavior. In 
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Figure 6. Plane wave excitation and re-radiation of the wave 
guided by a dielectric layer with periodically 
modulated dielectric constant for the case: a) plane 
wave incident from the right; and b) plane wave 
incident from the left. 
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Section IV, it is shown that total reflection can be achieved, and that 
the frequencies of total reflection can be predicted from the properties 
of the leaky wave inodes. 
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Ill . Formulation of the Reflection Problem 


In this section we develop the mathematical formalism for the re- 
flection and transmission coefficients at a layer of a periodic dielec- 
tric in a way that can be implemented on a computer for numerical evalua- 
tion. The analysis is restricted to the case of TE polarization (E x = E 2 
= 0) but it can be readily extended to the TM polarization. Expressions 
for the fields in the infinite periodic medium of Figure 2 are first 
derived. Subsequently the scattering matrix for a single surface normal 
to z, as shown in Figure 3, is found. Finally, network concepts are used 
to join the scattering matrices for the two surfaces normal to z that are 
shown in Figure 4, so as to obtain the reflection and transmission coef- 
ficients for the layer. 

A. Fields in an Infinite Periodic Medium 

In studying the fields in the infinite periodic medium we use the 
approach of Collin [11] and Lewis and Hessel [12] which expresses the 
fields as a superposition of modes having different wavenumbers along z. 
Assuming a time dependence exp(-ia»t), the fields of the n-th mode have 
dependence along z given by exp(+i* n z) . Thus the only non-zero component 
of electric field Ey and the z-component of magnetic field H z can be 
written as the sum of modal fields in the form 


E y (x,z) = Z [A n exp(i/e n z) + B n exp(-i* n z)] e n (x) , 
n=-« 

(3) 

Oft 

H z (x,z) = £ [A n exp(i/c n z) - B n exp(-i/c n z) ] h n (x) . 

n=-« 

Here e n (x) and h n (x) are the x-dependent mode functions while A n and B n 
are the modal amplitudes for waves propagating in the +z and -z direc- 
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tions, respectively. 

Because each slab of dielectric has a homogeneous e, the mode func- 
tions can be expressed in trigonometric form. We first define the wave- 
numbers along z in the two dielectrics as 


u n 88 < k l 2 ~ *n 2 >*' 
v n = (k 2 2 - « n 2 )*, 


(4) 


where kj 2 = k Q 2 s 1 and k 2 2 = k 0 2 e 2 . We further define the impedances and 
admittances for the dielectrics by 


z ln 88 1 ^ Y ln 88 ®^o/ u n' 
z 2n 88 l/*2n 88 


(5) 


Referring to the coordinate system in Figure 2, the mode functions 
for -dj , < z < 0 are given by 


e n (x) = V n cos^x) + i Z ln I n sin(u n x), 
h n< x > " J n cos(u n x) + i Y ln V n sin(u n x) . 

In the range 0 < z < d 2 , the mode functions are 


( 6 ) 


e n (x) * V n cos(v n x) + i Z 2n I n sin(v n x), 
h n (x) ■ I n cos(v n x) + 1 Y 2n v n sin(v n x) . 


(7) 


The constants V n and I n in (6) and (7) are defined using the Floquet 
condition discussed below. 

The Floquet condition requires that the fields at one end of a 
period (z = -dj ) differ from those at the other end (z = d 2 ) by at most a 
phase factor. Writing the phase factor as exp(iSgd), the Floquet condi- 
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tion is 


e n (d 2 ) • e n (-di) exp(iS 0 d), 

( 8 ) 

h n (d 2 ) = h n( -d l) exp(iS 0 d). 

Substituting from (5), (6) and (7) into (8) gives two homogeneous equa- 

tions in the two unknowns V n and I n . In order to have a non- trivial 
solution of these equations, it is necessary for u n and v n to satisfy the 
secular equation 


cos S 0 d = costUj^) cos(v n d 2 ) 

+ E( u n /V n> + < v n /u n>} sin(u n d 1 ) sin(v n d 2 ). (9) 

Since Ujj and v n are functions of « n , (9) can be viewed as a relation 
between « n and S 0 . As will be seen later, S 0 * k Q sinfl where 9 is the 
angle of incidence in Figure 1. Thus (9) serves as an equation whose 
roots are the allowed values of * n , and represents the dispersion equa- 
tion of the periodic medium. For large |n|, the roots are well ap- 
proximated by 

« n * [k 0 2 e a - (S 0 +n2«/d) 2 ]*, (10) 

which holds even for relatively small values of |n| . From expression 
(10), it is seen that only for small values of |n| will « n be real, 
whereas higher-order solutions will be below cutoff along z. 

When the dispersion equation (9) is satisfied, the ratio I n /V n can 
be determined form (8). After some manipulation, it is found that 


I n cos(v n d 2 ) - exp(iS 0 d) cos(u n d 1 ) 

= 1 Z 2n sin(v n d 2 ) + Z ln exp(iS Q d) sin(u n d 1 )' (11) 

whereas I n can be found from (11) if V n is known. The value of V n is 
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itself arbitrary and is usually obtained by normalizing the mode func 
tions. In this analysis, we use the normalization 


do 

/ I e ( x ) I l 2 dx = d . 
-*1 


(12) 


This normalization is carried out numerically during the computations, as 
described subsequently. 

Because of the Floquet condition, the mode functions e n (x) and h n (x) 
are periodic functions of x multiplied by the phase factor exp(iS Q x). 

Thus we may write e n (x) as the Fourier sum 


where 


e n (x) = l a nq exp(iS q x), 

qs-flp 

Sq “ Sq + q2s/d. 


(13) 


(14) 


Alternatively, the expansion coefficients a nq can be found from the 
integral 


l nq = d f e n (x) exp(-iS q x) dx. 

-di 


(15) 


Substituting (6) and (7) into (15), we obtain after much manipula- 
tion that 


n q 


(1) 
v n t J nq 


( 2 ) 

' J nq ) - 


(16) 


where 

( 1 ) 

J nq * i ( (S q + Y n ) [exp(iS q d 1 JcostUjjd! ) - 1] 

+ [u n + (Y n S q /u n )]exp(iS q d 1 )sin(u n d 1 ) } /(u n 2 - S q 2 ), 

with 
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iu n v n tcos(v n d 2 ) - exp( iS 0 d) cos ) ] 
u n sin(v n d 2 ) + v n ex P( iS Q d) sin(\i n d 1 ) 


(18) 


Y 


n 


( 2 ) ( 1 ) 

The quantity J ng In (16) is of the same form as that of J ng with u n 

replaced by v n , and dj replaced by -d 2 in (17). 

The normalization (12) is equivalent to requiring 


E I 


q=-co 


*nq' 


2 


1 . 


(19) 


Prom (16) it is seen that the normalization condition (19) implies 


r n 


( 1 ) ( 2 ) 
'nq “ J nq 


{ 2 I *^na “ ^na I ^ ^ • 

qa-< 


( 20 ) 


Substituting expression (13) into (3) for E y (x,z) and changing the 
order of summation gives 


E y (x,z) • I I [A n exp(i* n z) + B n exp(-i* n z)] a gn exp(iS_x), (21) 

qa-oo H M 


When applying boundary conditions at the surface z=0 in Figure 2, it 
is necessary to consider the x component of magnetic intensity H x (x,z). 
For the TE polarization, and using the Maxwell curl equations, it is 
readily seen that H x can be found from the derivative with respect to z 
of Ey. The resulting expression is 

-i®p 0 H x (x.z) 

CO oo 

58 E E * n t A n ®xp( 1 * n 2 ) “ B n exp(-i* n z)] a qn exp(iS x). (22) 

qa-oo n=-oo 
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Expressions (21) and (22) are used in the next section to find the scat- 
tering matrix for a single interface z ■ 0. 

B. Scattering at a Single Interface 

A single interface at z=0 is depicted in Figure 3. The field in the 
air region z < 0 consists of an incident plane wave propagating at an 
angle 0 with respect to the z axis, and reflected plane waves correspond- 
ing to the specular and higher space harmonics. We define 


S Q = k Q sin 0, 
C Q =« k 0 cos $ , 


(23) 


and assume the incident electric field to be polarized along y with 
amplitude E Q . The electric field of the incident wave is then given by 

E 0 exp (iS Q x) exp(iC Q z). (24) 


The electric field due to the reflected wave is the sum of the space 
harmonics and takes the form 

oo 

E R q exp(iSqX) exp(-iC q z) , (25) 

q=-oo 


where Sq is given by (14) and 

c q * (k Q 2 - S q 2 )*. (26) 

For the conditions of interest here S q 2 > k Q 2 for all q ^ 0 so that only 
the fields of the specular (q = 0) space harmonic propagate away from the 
interface, while the fields of all other space harmonics decay. The 
amplitude coefficients R q have yet to be determined. 

From (24) and (25) the total electric field in the air region is 
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seen to be 


E y (x,z) = E [<5 q o E o exp(iC Q z) + R q exp(-iC q z)] exp(iS q x), (27) 

qa-co 

where <5 q0 is the Kronecker delta. The x component of magnetic intensity 
in the air can be found from the derivative with respect to z of (27), 
which yields 

-i<ap 0 H x (x,y) 

* l c q [$ q0 E q exp(iC Q z) - R q exp(-iC q z)] exp(iS q x). (28) 

q=— oo 

The boundary conditions at z-0 require that E y and H x be continuous 

there. Equating the pair (21), (27) and the pair (22), (28), and making 

use of the orthogonality of the functions exp (iS q x) over a period, one 
obtains 

oo 

5 qO E o + R q " 2 [A n + B n ] a qn , (29) 

n 3 -» 

00 

c q <Wo - Rq) - l «n ' A n ‘ B nl a qn- « 30 > 

n s -« 

In (29) and (30), n ranges over all positive and negative integers, so 
that these equations represent two infinite sets of equations. 

We wish to solve (29), (30) for the amplitudes R q and A n of the 
waves traveling away from the surface in terms of the amplitudes E Q and 
B n of the incident waves. In this way, we obtain the scattering matrix 
of the surface. The waves incident from the periodic medium arise from 
reflection at the second surface, as suggested in Figure 4. For the 
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conditions of interest, only one or two waves are propagating in the 
periodic medium, while all other waves are cutoff along z. It can be 
shown that the two propagating waves correspond to the indices n = -1, 0. 
Since the two surfaces are separated by at least one half -wavelength, the 
fields of the cutoff waves excited at one surface are exponentially small 
at the other surface. Hence, to a good approximation we may assume B n = 

0 for n j* -1, 0. 

While all the higher modes are excited in the air and in the pe- 
riodic medium, the interaction at the two surfaces and radiation into the 
air are described by the amplitudes R Q , A Q , A_j^ of the propagating waves. 
Thus we ultimately need only the 3x3 portion of the full scattering 
matrix relating R Q , A 0 , A_ x to E 0 , B 0 , B_ x . This scattering relation 
takes the form 


o 


s i,i 

S l,0 S l,-1 


r 

E o 

A o 

A -i 

at 

s o , 1 
s -i,i 

S 0,0 s 0,-l 

S -l,0 S -1 , -1 

- 


®0 

B-l 


(31) 


To solve for the elements Saj5 in the scattering relation (31), we 
first multiply (29) by C g and then add (29) and (30). Since B n = 0 for 
n 9 * -1,0, the resulting equation may be written in the form 


2C 0 5 q0 E o + ( *0 “ c q) a 0q B 0 + (*-l " c q) a -lq B -l 

CO 

■ 2 (*n + c q> a qn A n* (32) 

n*-« 

Expression (32) represents an infinite set of equations with index q in 
an infinite number of unknowns A n . 

The terms ( * n + C g ) a gn can be viewed as elements of a matrix. With 
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this view, it can be shown that the elements decrease as one moves away 
from the main diagonal. Thus it is reasonable to solve (32) by using a 
finite truncation of the summation and a corresponding limitation on the 
number of values of g considered. For the parameters chosen in this 
study, it was found to be sufficient to allow q, n to range over the 
integers -3, -2, -1, 0, + 1, + 2. With this truncation, (32) is solved 
for A n with n* -3 ... +2. Returning to (29), we then compute R 0 from 

2 

R 0 3 " ®o + £ A n a Qn + a 00 B 0 + a o,-i ®-i* (33) 

n»-3 

Collecting the resulting values of A 0 , A_^ and R 0 due separately to E Q , 

B q and B_^ gives the values of the scattering matrix in (31). 

C . Scattering From A Periodic Laver 

Having found the scattering matrix for a single surface (z * 0), 
network concepts can be used to treat the interation with a second sur- 
face at z * h. As discussed previously, the interations between the two 
surfaces are essentially due to the propagating modes in the periodic 
medium. For our case, these are the n » -1, 0 modes, which are shown in 
the transmission line model for the interation shown in Figure 7. At the 
surface z * 0 the incident waves E Q , Bq and B.j couple to the scattered 
waves R 0 , A 0 and A.j . 

At the surface z * h, the incident waves in the layer are 
A 0 exp(i/c 0 h) and exp(i/c_ 1 h), while no wave is incident from the air 

side. In this case, the scattered waves are T 0 , B 0 exp(-ijc Q h) and 
B.j exp(-i*_ 1 h). The relation between scattered and incident waves is 
again given by (31), which for the foregoing conditions takes the form 
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Figure 7. Microwave equivalent network for determining reflection 
transmission at a periodically varying dielectric layer. 



( 34 ) 


T 0 


s l,l s l,0 s l,-l 


0 

B 0 exp(-iK 0 h) 

= 

S 0,l s 0,0 S 0,-l 


A 0 exp(i/t 0 h) 

B_ 1 exp( -i/c.jh) 


_ S -1 , 1 s -l,0 S -l,-l_ 


A_2exp( i/c^h) 


For the problem of scattering by a layer, the field E Q is known and 
one wishes to solve for the reflected and transmitted wave amplitudes R 0 
and T 0 , respectively. To this end, (31) and (34) can be viewed as six 
inhomogeneous equations in six unknowns Rq, Aq, A.j, Bq, B_j, Tq . Assum- 
ing E Q =1 , the solution of these equations for Rq and Tq give the reflec- 
tion and transmission coefficients of the layer. This approach has been 
used as the final stage in our computer program, as discussed below. 

D. Computer Program for R 0 and T 0 

A program has been written in the PL-1 language to compute Rq and T 0 
by the methods derived above. The listing of the program is given in 
Appendix A. An outline of the program is given below. 

1. Given input frequency at, angle of incidence 9, geometric parameters 
of the layer d 1# d 2 , h and the electrical parameters and e 2 . 

2. For integers n between -3 and 2, compute * n from (9) using Newton's 
method with starting value given by (10). 

3. For integers n, q between -3 and 2 compute A n g using (16) -(18) and 
( 20 ) . 

4. Using (32) and (33), solve for A n (-3 < n <2) and R 0 for E Q = 1 , Bq 
= B _2 = 0. This gives the elements Sj^ S Q j, S_ lfl in (31). 

Repeat for Bq = 1 with E Q = B.j = 0, and then for B_ 1 = 1 with E Q = 
Bq = 0 to get the remaining scattering coefficients in (31). 

5. Using (31) and (34) with E Q = 1, solve for Rq and Tq. 
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Several checks were carried out to ensure that the program was 
working properly. Choosing and e 2 very close to each other, we com- 
puted Rq and T q. As expected, they were very close to the values for a 
homogeneous dielectric layer of value e a . For values of and e 2 used 
subsequently, it was found that the scattering matrix in (31) conserves 
power. Finally it was observed that I Rq I 2 +IT Q I 2 is very close to unity, 
as required by power conservation. 
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IV . Numerical Studies of Frequency Selective Reflection 

Using the computer program described previously, we have carried out 
numerical studies for several examples. The purpose of these studies is 
to gain insight into the frequency selective behavior that can be ex- 
pected for Rq and Tq, and to relate this behavior to wave processes in 
the periodic layer. We have therefore arbitrarily chosen the angle of 
incidence 9 = 45° and set = d 2 ** d/2. 

For the initial studies, we have assumed ej = 2.56 and e 2 = 1.44, 
which are realistic values for low-loss plastics. With these choices, 
the average dielectric constant is e a ■ 2. The perodicity d must now be 
chosen such that, at the high frequency of interest, two modes (with n * 
-1, 0) propagate- in the dielectric layers. Furthermore, only the spe- 
cular (q = 0) space harmonic must propagate in the air. 

A. Choice of d, h and frequency 

The restriction on d needed to insure that only the q = 0 space har- 
monic propagates in the air can easily be interpreted with the help of 
Figure 8a. The incident wave has wavenumber S Q = k Q sin 9 < k Q along x. 
Wavenumbers S q of other space harmonics lie at a distance q(2n/d) away 
from S Q , as shown for q = -2, -1 and + 1 in Figure 8a. Provided that d 
is small enough so that 

k Q sin 9 - 2n/d < - k Q , (35) 

the q * -1 space harmonic will lie outside the visible circle defined by 
S 2 + C 2 = k Q 2 . It is further seen that, if (35) holds, all other space 
harmonics will also lie outside the circle, so that C q defined by (26) is 
imaginary for q ? 0, and the space harmonics decay away from the layer. 
Provided that the modulation (Sj - s 2 )/e a is not too large, the concept 
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(a) 



Figure 8. Vissible circle for determining the propagating space harmonics 

in: a) air; and b) a medium with dielectric constant e . 

a 
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of the visible circle can also be used to estimate the number of propa- 
gating waves in the layer. When the modulation is small, (10) can be 
used as an estimate for « n , in which case the visible circle is given by 
S n 2 + /c n 2 = k 0 2 e a , as shown in Figure 8b. It is seen from this figure 
that, for the two modes n* -1, 0 to propagate in the layer, d must satis- 
fy 


k Q sin 9 - 4n/d < -k Q /e^ < k Q sin 9 - 2n/d, 
k Q /e^ < k Q sin 9 + 2?c/d. 


(36) 


Conditions (35) and (36) can be rearranged into the following in- 
equalities : 


d/A < 1/(1 + sin 9 ) , 

+ sin 9) < d/A < 2/(JT ~ + sin 9), (37) 

d/A < l/(/!Tj[ - sin 9). 

Assuming e a ■ 2 and 9 = 45°, these inequalities are d/A < 0.586, 

0.471 < d/A < 0.943 and d/A < 1.414, respectively. To satisfy these in- 
equalities, we have chosen d/A - 0.54. 

Initially a value of layer thickness h at which total reflection 
will occur was obtained by computing Rq for various values of h /A. 
Subsequently, it was found that values of h and d for total reflection 
could be related via the conditions for guidance of a wave by a layer of 
uniform dielectric constant e a . Whereas our initial approach gave us the 
value of h / A = 0.925 for sample calculations, it is the subsequent inter- 
pretation that is discussed below. 
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B. Variation of R 0 With Frequency 


Computations of the frequency dependence of Rq have been made assum- 
ing that h = 0.925 and d = 0.54 for e-^ = 2.56 and e 2 = 1.44. This choice 
produces total reflection for a frequency f such that X ~ c/f is about 
unity. Note that, if h and d are scaled by X, then total reflection can 
be obtained at any desired frequency. The results of the calculation for 
I Rq I are depicted in Figure 9, where we have used the normalized frequen- 
cy variable k Q h = 27tfh/c, and have plotted up to the value k Q h = 6.30 at 
which the q = -1 space harmonic in air switches from cutoff to propagat- 
ing along z. 

For k Q h < 5.12, only the n = 0 mode in the periodic dielectric is 
propagating along z. In the frequency range 0 < k Q h < 5.12 for single 
mode propagation, R 0 vanishes at the two frequencies k Q h = 2.56 and 5.03, 
at which « 0 h * 3.145 and 6.355. Thus, frequencies of total transmission 
occur when the layer thickness is close to a multiple of one half the 
effective wavelength along z, as predicted in Section II. The difference 
between the values of /c Q h for total transmission and n , 2k are due to the 
non-zero phase of the transmission and reflection coefficients at the 
individual surfaces z * 0 and h, which result from excitation of higher 
cutoff space harmonics. 

For k Q h < 5.12, where the n = -1 mode also propagates along z in the 
layer, total reflection takes place at two frequencies given by k Q h = 

5.32 and 5.83. In the vicinity of these frequencies for total reflec- 
tion, the variation of I Rq I is that associated with resonances wherein a 
frequency dependent function has a real axis zero and a nearby pole at a 
complex location. 
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To examine the bandwidth of the total reflection resonance, we have 
plotted I Rq 1 2 on an expanded scale in Figure 10. Since IR 0 I 2 + ! Tq I 2 = 

1, the curve can also be used to determine I Tq 1 2 by using the vertical 
scale to the right of the plot. The wider of the two peaks of I Rq 1 2 is 
centered at k Q h = 5.83. For this peak, the fractional bandwidth between 
the frequencies at which I Rq 1 2 * 0.9 is 0.8635. 

The narrower of the two peaks in Figure 10 is shown further expanded 
in the insert. For this peak, the fractional bandwidth between the fre- 
quencies at which I Rq 1 2 = 0.9 is less than 0.0435. By comparison, much 
wider bandwidths are found for total transmission when the n = -1 mode in 
the layer is cut off. For example, in a region about k Q d * 2.56 in 
Figure 9, a 4035 bandwidth is found between the frequencies for which 
I R 0 1 2 « 0.1. 

C . Prediction bv Means of Guided Waves 

The location of the frequencies of total reflection can be predicted 
from the properties of the waves guided by the periodic layer. Consider 
first the case of a wave guided along a uniform layer having dielectric 
constant e a . The normalized propagation constant 0gh of this guided wave 
is plotted horizontally in Figure 5, versus the normalized frequency k Q h, 
which is plotted vertically. The lowest g = 0 guided wave mode starts at 
the origin and becomes asymptotic to the wavenumber k 0 /s^ of the layer. 
Higher guided-wave modes start at points k Q h = gjr(e a -l) along the 45° 
line, where g = 1, 2, 3 .... 

In Figure 11, we have repeated the plot of Figure 5, and have added 
the dispersion curves for waves propagating in the negative x direction 
(8 < 0). We have also plotted as a broken line the transverse wavenumber 
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Dispersion curves for surface waves traveling in the ±x direction (solid curves) and 
the normalized transverse wavenumber S h of the incident plane wave (broken line). 




S Q h = k Q h sin 9 of the incident plane wave. Except at k Q h = 0, the 
broken line is never close to the dispersion curves for the guided waves. 
As a consequence, an incident plane wave cannot couple to the guided 
waves on a uniform layer. 

If the dielectric layer is made periodic along x, the field of each 
guided wave mode becomes a sum of space harmonics. For the guided waves 
traveling in the -x direction, the wavenumber of the q = -1 space har- 
monic is (~0g + 2*/d) . When normalized by h, the dispersion curve of 
this space harmonic has the same form as -0gh versus k Q h, except for a 
shift 2nh/d to the right. While finite modulation affects the value of 
0g for the guided wave, for small modulation of the dielectric constant 
0g is close to that for a uniform layer. 

Dispersion curves for the q = -1 space harmonics of the guided waves 
in the small modulation limit are shown in Figure 12 for h/d = 

0.925/0.54. We have also drawn a broken line representing S Q h versus 
k Q h. Intersection of the S Q h line with the dispersion curvers indicates 
strong coupling between an incident plane wave and guided waves through 
the q = -1 space harmonic. Note that, above the dashed line having an 
angle of -45°, the q » -1 space harmonic in the air propagates along z. 
Thus for reflection and transmission of a single space harmonic, the 
operating point along the S Q h line must be kept below the dashed line. 

For the parameters used in drawing Figure 12, this condition implies that 
k Q h < 6.30 for one propagating space harmonic in air. 

In Figure 12, the line S Q h intersects the dispersion curve for the g 
= 0 guided wave at k Q h * 5.27, and for the g = 1 guided wave at k Q h = 
5.72. These values are close to the values k Q h * 5.32 and 5.83 for total 
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Figure 12. Transverse wavenumber S Q h of the incident plane wave and the dispersion 
curves of the q=-l space harmonics for the g=0,l,2 guided waves in the 
limit of small modulation for e =2 and h/d=0. 925/0.54. 



reflection ( 1 Rq 1=1) obtained from Figure 11. The deviation between the 
values of k Q h obtained from Figures 11 and 12 is thought to result from 
the fact that the finite modulation of the dielectric constant of the 
layer alters /9g from the value obtained for a uniform slab. Thus, de- 
creasing modulation should bring the values closer together, while in- 
creasing modulation should result in greater deviation. This latter 
condition is shown subsequently. 

To further demonstrate the relation between the frequencies of total 
reflection and the guided waves of the layer, we have considered a layer 
of increased thickness h = 1.1, but the same periodicity d = 0.54. The 
dispersion curves of the q * -1 space harmonics of the first three guid- 
ed-wave modes are shown in Figure 13 for the limiting case of small 
modulation. The broken line giving S Q h = k Q h sin 0 is seen to intersect 
the three dispersion curves at k Q h * 6.22, 6.67 and 7.30. Our model pre- 
dicts that total reflection should take place at normalized frequencies 
close to these values. 

A plot of I R 0 I versus k Q h for h = 1.1 and d = 0.54 is shown in 
Figure 14. From this plot, total reflection is seen to occur at the 
three frequencies k Q h =6.25, 6.78 and 7.42, which are close to those 
predicted by the small modulation theory. The bandwidth over which I Rq I ^ 
>0.9 about each frequency of total reflection is seen to increase as k Q h 
approaches the value 7.50 where the line S Q h crosses the dashed line, 
above which the n = -1 space harmonic propagates in air. The bandwidth 
about the lowest of these frequencies is only 0.01%, while that of the 
highest is 0.7%. 
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Figure 13. Transverse wavenumber S o h of the incident plane wave and the dispersion 
curves of the q=-l space harmonics for the g=0,l,2 guided waves in the 
limit of small modulation for e =2 and h/d=l.l/0.54. 




0*1 
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Figure 14. Variation of |Rq| with normalized frequency k Q h for e^=2.56, £ 2 =! . 44 and h/d=l.l/0.54. 



c. 


Influence of Modulation 


To explore the influence of modulation, we have computed the reflec- 
tion and transmission coefficients for a layer with ©2=3 and e 2 = 1- 
This layer has average dielectric constant e a = 2, as before. We further 
assume that d = 0.54 and h =0.925, as in the case of the results pre- 
sented in Figures 10-12. A plot of IR Q I versus normalized frequency k Q h 
is shown in Figure 15. The variation of IRqI is seen to be qualitatively 
the same as that of Figure 10. The increased modulation is seen to shift 
the first frequency of total reflection (IRqI = 1) to k Q h = 5.45 and the 
second to k Q h = 6.12, which are farther from the respective values 5.27 
and 5.72 predicted by small modulation theory. 

Besides shifting the frequency of total reflection, the modulation 
influences the bandwidth. At the first total-reflection frequency, the 
bandwidth for IRqI 2 > 0.9 is 0.001*. However, at the higher total re- 
flection frequency, the bandwidth is 2.5*. The modulation is also seen 
to have a small effect on the frequencies of total transmission (IRqI=0). 
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V. 


Conclusion 


It has been shown that frequency selective reflection and transmis- 
sion takes place at a periodically modulated dielectric layer. Frequen- 
cies of total transmission and total reflection were found, and they can 
be related to various wave phenomena. In the limit of small modulation, 
these frequencies can be estimated from the appropriate wave phenomenon 
in a uniform layer having dielectric constant equal to the average of 
that in the periodic layer. 

In the range of low frequencies where a single space harmonic propa- 
gates along z in the periodic dielectric, total transmission occurs when 
the layer thickness h is one half the effective wavelength along z, i.e., 
when h * k/kq. For small modulation, kq » ( e a k 0 2 -S 0 2 ) **, where e a is the 
average dielectric constant. Total reflection can be achieved at those 
higher frequencies for which two space harmonics propagate along z in the 
periodic dielectric. These frequencies of total transmission are asso- 
ciated with the excitation of leaky waves guided by the dielectric layer. 
In the limit of small modulation, the frequency of total reflection can 
be approximated from the dispersion characteristics of waves guided by a 
uniform dielectric layer. 

In the examples treated, the bandwidth over which IRqI 2 > 0.9 about 
the frequency of total reflection was found to be small. The largest 
bandwidth obtained was 2.5%. While angle sensitivity was not computed, 
the narrow frequency bandwidth suggests that, at the frequency of total 
reflection, I Rq I 2 will be sensitive to the angle of incidence e. 

Whereas the study was carried out only for the TE polarization, the 
form of the results have implications for the TM polarization. We expect 
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that the frequencies of total reflection for the TM polarization are also 
associated with the excitation of the leaky waves guided by the periodic 
layer. However, the dispersion characteristics of the leaky TM waves 
will differ from those of the TE polarization. As a result, it is ex- 
pected that incident plane waves of the TE and TM polarizations will, in 
general, experience total reflection at different frequencies. Hence, 
the periodic dielectric layer is expected to be polarization sensitive. 
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//TEMNNPOF' JOB (F 109602, CHEO) , 'L.S.CHEO' 

// EXEC PLIXCLG 
//PL1L.SVSIN DD * 

TEMODE : PROCEDURE OPT I ONS (MAIN); 

/* FINAL PROGRAM FOR TE MODE WITH ONE AND TWO */ 

/* PROPAGATION MODES */ 

/* GAMA(-l) CAN BE EITHER REAL OR IMAGINARY */ 

/********************-*#-****************■***■****** *■****•*/ 

/* INSTRUCTION FOR DATA ENTRY : */ 

/* THE FIRST DADA CARD ENTRY ORDER IS : */ 

/* D,D1 ,D2, El ,E2,LAMTA, DEGREE, H */ 

/* THE SECOND DATA CARD IS FOR VALUES OF M AS : */ 

/* Ms-6,-5,-4,-3,-2,-1,0,1,2,3,4,5 */ 

/************-*******************-**********************/ 

/* CONSIDER V (-6: 5) ,U (-6: 5) , 5N(-6:5) AS ARRAYS */ 

/* XI = U/V, U IS REAL , V CAN BE REAL OR IMAGINARY */ 

/* */ 

/* FIND A ROOT NEAR UO=SO OF THE FUNCTION */ 

/* F (U) = 0 BY NEWTON'S METHOD , TO SIX DIGIT */ 

/* ACCURACY. */ 

/* */ 

/* THIS PROGRAM COMBINES TEMNNPO , RTLCOMP , AND RTLGAUS */ 

DCL < AA (1:8,1: 11) ,X(1:3,1:8) ) FLOAT ; 

/* AA IS THE COEFFICIENT MATRIX FOR GAUSS ELIMINATION */ 

/* AND X IS THE SOLUTION VECTOR OBTAINED */ 

/* FROM THE GAUSS ELIMINATION METHOD PROCEDURE */ 

/* */ 

/* MATRIX MN ( M , N > , EXCEPT PRINTING PART FOR ALPHA*/ 

/* SUMJREAL , SUMJ IMAG , SUMJ1J2 ARE ELIMINATED. **/ 

/*** ***/ 

DCL <S,SS,U, (U1,U2) (-6:5) , EPSILON) FLOAT DEC; 

DCL ( V , SM) (—6: 5) FLOAT , 

(K , KMAX ) FIXED DEC; 

DCL (UD1 (—6: 5) , RTRD2 , D2 , D1 , D , SO , SQRTR) FLOAT ; 


DCL (El , E2,K0 , SOD , UR , RU , SUM1 , SUM2 , R (-6; 5) ) FLOAT ; 

DCL (A, B, TERM, PI ,SMD(-6:5) , GAMAIN , GAMA (-6: 5) ) FLOAT ; 

DCL M FIXED (4 , O) INITIAL (-5) ; 

DCL (Y,YR,YI ,Y1 , Y2 , Y3 , Y4 , YMAG) (-6:5, -6:5) FLOAT ; 

DCL (J1R1,J1R2, J1R3,J1R4, J1R) (-6:5, -6:5) FLOAT; 

DCL <J1I1,J1I2,J1I3, J1I4,J1I) (-6:5, -6:5) FLOAT; 

DCL (J1 , J2) (-6: 5,-6: 5) FLOAT; 

DCL ( J2R1 , J2R2 , J2R3 , J2R4 , J2R) (-6: 5, -6:5) FLOAT; 

DCL (J2I1 ,J2I2,J2I3,J2I4,J2I) (—6: 5 , —6; 5) FLOAT; 

DCL (SUMJ 1J2,SUMJR, SUMJ I ) (-6; 5, —6: 5) FLOAT; 

DCL < AMN , AMNR , AMN I , MN , MNR , MN I ) (-6:5,-6:5) FLOAT; 

DCL SUML (—6: 5) FLOAT ; 

DCL (ALPHA, CN,C0SVD2, SI NVD2) (-6:5) FLOAT; 

DCL (UMD1 , SND , SND 1 , SND2 , VMD2) (—6:5) FLOAT; 

DCL (C0SHVD2,SINHVD2) (-6:5) FLOAT; 

DCL <TR,TI) (3,-2; 1) FLOAT (6); 

DCL ( SCATR , SCAT I ) (-2: 1,-2: 1) FLOAT (6); 

DCL (ROR,ROI) (1:3) FLOAT (6) ; 

DCL ( L AMTA , DEGREE , THETA , K02 , RAT I 0 ) FLOAT ; 

DCL ( N , I , J , P ) F I XED (5,0) ; 

DCL ( RTLR (4,4) , RTL I (4,4) ) FLOAT; 

DCL ( AB ( 12) , CD ( 12) ) FLOAT ; 

DCL GAMA0_1H FLOAT ; /* GAMA (0) *GAMA (-1 ) *H */ 

DCL GAMA02H FLOAT; /* 2*GAMA(0)*H */ 

DCL GAMA 12H FLOAT: /* 2*GAMA(-1)*H */ 
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DCL GAMAOH FLOAT ; /* GAMA (0) *H */ 

DCL GAMA_1H FLOAT; /* GAMA(-1)*H */ 

DCL (LL) FIXED (5,0) ; 

DCL FLAG2 FIXED (2,0); 

DCL ( RREAL , R I MAG , RPR I MER , RPR I ME I ) FLOAT ; 

DCL (CDISC , GAMAD) (-6:5) FLOAT; 

DCL ( B2R , B2I , B3R , B3I ) (-6:5) FLOAT;' 

DCL H FLOAT ; 

DCL (MM , NN , PIVOT , RS) FIXED (5,0) ; 

DCL ( TOR , TO I , TM I NUS 1 R , TM I NU5 1 1 > FLOAT ; 
DCL ( RABS , RPHASE , RPRABS , RPRPH ASE ) FLOAT ; 


/* MM=NO . OF COLUMNS & NN=NO. OF ROWS IN MN (N , M) •*/ 

/* MATRIX IN PROGRAM TEMNNPO */ 

/* IF PIV0T=1 THEN THE PIVOT IS SUBSCRIBED, */ 

/* RS=NO. OF RHS OF THE AUGMENTED COEFFICIENT MATRIX*/ 

/* IN THE PROCEDURE GAUSS (ELIMINATION METHOD */ 

/* */ 

/* DEFINE THE FUNCTION F(X) */ 


F ; PROCEDURE (U,M) 

DCL U FLOAT ; 

DCL M FIXED (4,0) ; 

CALL CHECKR (U , M) ; 

S^COS (UD1 (M) ) *A— 0 . 5*SUM 1 *S I N ( UD 1 (M) ) *B-COS (SOD) ; 

RETURN (S) ; 

END F; 

/■* */ 
/* DEFINE THE 1ST DERIVATIVE OF F(X> */ 

FPR I ME : PROCEDURE (U,M) ; 

DCL U FLOAT ; 

DCL M FIXED (4,0) ; 

DCL ( TERM 1 , TERM2 , TERM3 ’> FLOAT; 

CALL CHECKR (U,M) ; 

TERM1 =HD1*SIN(UD1 <M) ) *A-D2*C0S (UD1 (M) ) *B*tJR ; 
TERM2=-0.5*SUM2*SIN(UD1 (M) ) *B ; 

TERM3=-0 . 5*SUM 1 * ( D 1 *COS ( UD 1 (M) ")*B+TERM> ; 

SS=TERM1 + TERM2 + TERM3 ; 

RETURN (SS) ; 

END FPR I ME; 

/* PROCEDURE TO CHECK THE SIGN OF R=U**2-K0**2* ( E 1 -E2 ) •*/ 

/* */ 
CHECKR s PROC (U, M) ; 

DCL U FLOAT ; 

DCL M FI XED (4,0) ; 

UD1 <M) = U * D1 ; 

R < M ) =U**2-K0**2* ( E 1 -E2 ) ; 

SQRTR » SORT ( ABS (R (M) ) ) ; 

RTRD2 = SQRTR*D2 ; 

RU =» SQRTR/U ; 

UR = U/SQRTR ; 

SOD = SO * D ; 

IF R (M) > 0 THEN 

DO ; 

A = COS (RTRD2) ; /* WHEN R > 0 */ 

B = SIN (RTRD2) ; 

SUM1 = RU + UR ; 

SUM2 = 2/SQRTR— RU/U-UR**2/SQRTR ; 

TERM = D2*SIN(UD1 (M) ) *A ; 

END ; /* END OF R>0 CASE */ 

ELSE 

DO ; /* WHEN R< O */ 

A - COSH (RTRD2 ) ; 

B = SINH (RTRD2) ; 

SUM1 = UR - RU ; 

SUM2 = 2/SQRTR+RU/U + UR**2/SQRTR ; 

TERM = ~D2*SIN (UD1 (M) ) *A ; 

END s /* END OF R<C CASE */ 
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END CHECKR ; /** END OF CHECKR PROCEDURE **/ 

/* */ 
FLAG=1 ; 

ON ENDFILE (SYSIN) FLAG-0; 

/* MAIN PROCEDURE */ 

PI « 3.14156 ; EPS I LON= 1 . OE— 06 ; KMAX-SO ; 

GET LIST (D , D1 , D2 , E 1 , E2 , LAMTA , DEGREE , H) ; 

KO = 2*PI /LAMTA ; 

PUT SKIP; 

IF LAMTA- 1.0 THEN 

PUT SKIP EDIT ( ' TE MODE WITH TWO PROPAGATION MODES') 

(X (5) , A) ; 

ELSE 
DO ; 

PUT SKIP ; 

PUT SKIP EDIT ( ' TE MODE WITH ONE PROPAGATION MODE') 

<X <5) , A) ; 

END ; 

PUT SKIP EDIT (REPEAT ( ' * ' , 40) ) ( X (5) , A) ; 

PUT SKIP EDIT ('D-',D,'D1=',D1 ? ' D2= ' , D2, ' THET A= ' , DEGREE) 

< X (5) ,4 (A,F<9,3> ,X (4) ) ) ; 

PUT SKIP EDIT ( ' E 1 = ' ,E1 , 'E2=' ,E2, ' LAMTA= ' , LAMTA, 'K0-' ,K0> 

< X (5) ,4 (A , F (8 , 5) , X (3) > ) ; 

PUT SKIP; 

PUT SKIP ED I T ( ' H= ' , H ) (X (5) , A , F ( 8 , 4 ) ) ; 

PUT SKIP EDIT (REPEAT ( ,50) ) (X (5) ,A) ; 

PUT SKIP ; 

CALL INPUT (M) ; 

LOOP: DO WHILE <FLAG=1 > ; 

CALL CALCULATE (M) ; 

CALL PRINT (M) ; 

FLAG2- 1 ; 

CHECK 1 : DO WHILE (ABS ( (U2 (M> -U1 (M) ) /U1 (M) ) >*EPSI LON & 
FLAG2= 1 > ; 

IF K<— KMAX THEN DO; 

U1(M)=U2(M>; 

CALL CALCULATE (M> ; 

CALL PRINT (M) ; 

END; 

ELSE DO; 

PUT SKIP (2) EDIT ( 'FAILS TO CONVERGE ')( A ) ; 
FLAG2-0; 

END; 

END CHECK 1; 

PUT SKIP (2) ; 

PUT SK I P ( 2 ) EDIT( 'R*' ,R(M> > (X(5) ,A,E(14,6) > ; 

IF R(M) < 0 THEN 

PUT SKIP (2) EDIT ( ' V IS IMAGINARY ')( X (5) , A) ; 

V (M) * SQRTR ; 

PUT SKIP (2) EDIT ( ' V (M) - ' , V (M) ) (X(5) , A f E ( 14 , 6) ) ; 

GAMA IN = K0**2*E1-U2(M>**2 ; 

PUT SK I P ( 2 ) ED I T ( ' K0**2*E 1 -U2**2= ' , GAMA IN) (X(5) ,A,E(14,6 
IF GAMA IN < 0 THEN 

PUT SKIP (2) EDIT ( 'GAMA IS IMAG I NARY ' ) ( X ( 5 ) , A ) ; 

GAMA (M) - SORT (ABS (GAMAIN) > ; 

PUT SK I P ( 2 ) ED I T ( ' GAMA= ' , GAMA (M) ) (X(5) , A , E ( 1 4 , 6 ) ) ; 

IF FLAG2=1 THEN CALL OUTPUT (M) ; 

CALL INPUT (M); 

END LOOP; 

/* */ 
/* TO CQMF : ’UTE UM*D 1 , VM*D2 , SN*D 1 , SN*D2 , SN*D **/ 

/* TO COMPUTE C0S(VMD2) ,C0SH(VMD2) ,SIN(VMD2) ,SINH(VMD2) **/ 
/**********************************#******.******.**..**.*.*.***.*/ 

LOOPM : DO M*-6 TO 5 ; 

UMD 1 (M) = U2(M) * Dls 
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VMD2 ( M) = V(M) * D2 ; 

SND 1 <M) = SM<M> * D 1 ; 

SND2 ( M) = SM(M) * D2 ; 

SND (M ) = SM<M) * D ; 

SINVD2(M> = SIN (VMD2 (M) > ; 

SINHVD2 (M) » SINH(VMD2(M> > ; 

C0SVD2 (M) = COS ( VMD2 (!i) ) ; 

C0SHVD2 (M) = COSH (VMD2 <M) ) ; 

END LOOPM ; 

/***************************************************^-«*^****/ 

/* TO COMPUTE Y(M,N> **4H**********************************/ 

/* */ 

LOOPYM : DO M = -6 TO 5 ; 

LOOPYN : DO N = -6 TO 5 ; 

Y1(M,N) * COS ( UMD 1 ( M ) ) *S I N < SND ( N > ) ; 

Y4(M,N) = SIN (UMD 1 ( M) ) * SIN (SND <N) ) /U2 (M) ; 

IF R (M) > 0 THEN 
DO ; 

Y3 ( M , N ) =S I NVD2 ( M ) / V (M) +SIN (UMD1 (M) )*COS<SND(N> ) /U2(M) 
Y2 (M , N) a =C0SVD2 (M) -COS ( UMD 1 (M) >*COS(SND(N) ) ; 

END; 

ELSE 
DO ; 

Y3 <M,N)=SINHVD2<M) /V <M> +SIN (UMD 1 (M) )*COS(SND(N) ) /U2(M) 
Y2 < M , N ) —C0SHVD2 ( M ) -COS ( UMD 1 (M) >*COS<SND(N) ) ; 

END ; 

/****^*********************************************-«-********/ 

’ /* TO COMPUTE ABS<Y(M,N>**2) **/ 

/* TO COMPUTE REAL AND IMAGINARY PART OF Y: YR , YI **/ 

YMAG (M , N) =Y3 (M , N) **2+Y4 (M , N) **2 ; 

YR (M , N> = ( Y 1 (M,N) *Y3<M,N) +Y4 (M,N) *Y2 (M,N) ) /YMAG (M , N) ; 

YI <M, N> = ( Y2 (M , N) *Y3 (M , N) — Y 1 <M ,N> *Y4 <M , N> ) /YMAG (M , N ) ; 

END LOOPYN ; 

END LOOPYM ; 

/******************************************■»('*******■«*■»<'******•*•*/ 

/* TO COMPUTE J1(M,N) */ 

/* REAL AND IMAGINARY PART OF JI s J1R, J1I */ 

PUT SKIP ; 

L00PJ1M s DO M = -6 TO 5 ; 

DO N = -6 TO 5 ; 

J1R1 (M,N)«YI <M,N)*( 1-COS (UMD 1 (M) >*C0S(SND1 (N) ) ) ; 

J1R2(M,N)«(SM<N)+YR(M,N) )*C0S(UMD1 (M) )*SIN(SND1 (N) ) ; 
J1R3<M,N)=(U2(M)+SM<N)*YR<M,N) /U2(M) >*C0S(SND1 (N) ) ; 

J 1 R4 ( M , N ) ®S I N < SND 1 (N) )*SM(N)*YI (M,N> /U2(M) ; 

J 1R (M , N) — ( J 1R1 (M , N) — J 1R2 (M,N) +SIN (UMD1 <M> )* 

( J 1R3 (M , N) — J 1R4 (M , N) ) ) / <U2 <M> **2-SM <N> **2) ; 

Jill (M,N) = <SM (N) +YR (M , N) )*(C0S(UMD1 (M) >*C0S(SND1 (N) )-l ) ; 
J1I2(M,N)*YI ( M , N ) *COS ( UMD 1 <M> )*SIN(SND1 (N) ) ; 
J1I3(M,N>=SM<N)*YI (M , N) *CQS (SND1 (N) ) /U2(M) ; 

J 1 1 4 ( M , N ) «S I N ( SND 1 (N) ) * ( U2 ( M ) +SM < N ) * YR ( M , N ) /U2<M) ) ; 

J1I (M,N)*(J1I1 (M , N) — J 1 12 (M, N) +SIN (UMD1 <M) >* 

(JlI3(M f N) + J 1 1 4 < M , N ) ) >/(U2<M)**2-SM(N>**2) ; 

END ; 

END L00PJ1M; 

/*************************************************** **********-■#•/ 
/* TO COMPUTE J2(M,N) */ 

/* REAL AND IMAGINARY PART OF J2 : J2R , J2I */ 

PUT SKIP; 

L00PJ2M : DO M * -5 TO 5 ; 

DO N = -j TO 5 ; 

J2R4(M,N)~SM(N)*YI (M , N ) *S I N < SND2 ( N ) ) /V(M) ; 

J2I3(M,N) «SM<N) *YI ( M , N ) *COS ( SND2 ( N ) ) /V(M) ; 

IF R (M) >0 THEN 
DO; 

J2R1 (M,N> =YI (M , N) * (C0SVD2 (M> *COS (SND2 (N) >-l) ; 

J2R2 (M . N) = (SM (N ) +YR < M . N > > *C0SVD2 (M) *-SIN (SND2 <N) ) : 
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J2R3(M,N) = (0(M)+SM(N)*YR(M,M)/V(M> > *C0S (SND2 (N) ) ; 
J2R(M,N) = < J2R1 < M , N ) - 0 2R2 ( M , N ) +SINVD2 <M> * 

(J2R3(M,N)+J2R4<M,N) ) > / <V <M> **2-SM (N> **2> ; 
02I1<M,N)=(SM<N>+YR<M,N> > * (C0SVD2 (M> *C0S (3ND2 <N) >-l) ? 
J2I2 (M , N) =C0SVD2 CM) *S IN ( SND2 ( N ) ) *YI (M,N> ; 
02I4(M,N)=SIN(SND2(N> > * ( V (M) +SM <N) *YR (M , N) /V <M) ) ; 

021 (M , N) = (-J2I 1 (M,N)-J2I2(M,N)+SINVD2<M>* 

(0213 (M , N) —0214 (M , N) ) ) / ( V (M) **2-SM (M) **2) ? 

END ; /*END OF R (M) > O FOR J2 COMPUTATIONS ***/ 

ELSE 
DO ; 

02R1 CM,N)=YI <M,N>* (C0SHVD2 <M> *CQS < SND2 (N> > -1 > ; 

J2R2 (M , N) = (SM (N ) +YR <M , N) ) *C0SHVD2 <M> *SIN (SND2 (N) > ; 
02R3(M,N>=(-V(M)+SM(N)*YR(M,N) /V(M) > *COS (SND2 <N) > ; 

02R (M , N) =— (02R1 <M ? N) — J2R2 (M , N) +SINHVD2(M>* 

( J2R3 (M , N) +J2R4 (M , N) ) ) / ( V (M) **2+SM (N) **2> ; 
J2I1(M,N)=(SM(N)+YR(M,N) > * (C0SHVD2 <M) *COS (SND2 <N) )-l) ; 
02I2(M,N)=C0SHVD2(M)*SIN(SND2(N> )*YI (M,N) ; 

0214 (M ,N) =SIN (SND2 (N) ) * (-V (M) +SM (N) *YR (M ,N> /V (M) ) ; 

021 (M,N)=-(-J21 1 <M,N)-J2I2(M,N> +SINHVD2 (M> * < J2I3 (M , N) 
—J2 14 (M , N) ) ) /<V(M)**2+SM(N)**2) ; 

END ; /* END OF R(M> <« 0 FOR J2 COMPUTATION ***/ 

END ? 

END LOOPJ2M ; 

/************************#**************#**********************/ 

/* */ 


/* TO COMPUTE SUM01J2 

/* REAL AND IMAGINARY PART OF SUMJ1J2 : SUMJR , SUMO I 
SUMJM s DO M = -5 TO 5 ; 

SUMJN : DO N = -5 TO 5 ; 

SUMOR (M, N) =J1R(M,N> +J2R(M,N> ; 

SUM J I <M f N)*JlI (M,N) +J2I (M , N) ; 

SUMO 1 02 (M , N> =SUMJR <M , N) **2+SUM0I (M,N) **2 ; 

END SUMON; 

END SUMJM ; 

/****** TO COMPUTE ALPHA (M), M=-6 TO 5 ***************/ 

/* TRUNCATE INFINITE SUM TO SUM SUMO 102 <M,N) */ 

/* */ 

ALFA : DO M - -5 TO 5 ; 

SUML (M> =0 : PUT SKIP : 

SUMLL : DO LL ■ -5 TO 5 ; 

SUML (M> = SUML (M) + SUMJ 1 02 (M , LL) ; 

END SUMLL? 

/***** TO COMPUTE ALPHA <M> ********/ 


/* 


*/ 


*/ 

•*/ 


ALPHA <M) « SORT (SUML (M> > ; 

PUT SKIP EDIT ( ' ALPHA ( ' T M , ' ) — ' , ALPHA (M) ) <X ( 5 ) ,A,F< 2 , 0 > V A,E< 12 , 5 > 


END ALFA; 

/** TO COMPUTE CN(N> , GAMA (M) , AND MN(N,M> M , N = -6 TO 5 **/ 

/** **/ 

PUT SKIP; 

PUT SKIP EDIT (REPEAT ( ,40) ) (X (3) ,A> ? 

PUT SKIP; 

LOOPCN; DO N * -5 TO 5; 

CDISC(N) = K0**2-SM(N>**2 ? 

CN (N) = SORT (ABS (CDISC (N) ) ) ; 

PUT SKIP EDIT ( # CN( ' ,N, ' > - ' ,CN(N) , ' CDISC ( ' ,N, ' > ,CD ISC (N) > 

( X (3) ,2 <A,F<3,0> ,A,E<12,5> ,X (2) ) ) ; 

END LOOPCN ; 

GAMADM: DO M = -5 TO 5 ; 

GAMAD(M) = K0**2*E1HJ2 (M) **2 ; 

GAMA <M)=SQRT (ABS (GAMAD <M) ) ) ; 

PUT SKIP EDIT ( ' GAMA ( ' ,M, ' ) = ' , GAMA (M ) , ' GAMAD ( ' ,M, ' >«' , GAMAD ( 
(X (3> ,2 (A,F (3.0) , A.E < 12.5) .X (2) ) ) : 
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OF PGuK QUALITY 

END GAMADM ; 

LOOPAMN: DO M = -5 TO 5; 

LOOPMNA : DO N = -5 TO 5* ; 

AMNR (M,N) -SUMJR ( M ,'N) / ALPHA (M) ; 

AMNI <M,N)=SUMJI <M,N> /ALPHA (M) ; 

IF CD ISC (N) >0 & GAMAD(M) > 6 THEN 

CASE 1 : DO ; 

MNR <N , M) = <CN<N)+GAMA(M> )*AMNR(M,N) ; 

MNI (N , M) = <CN<N)+GAMA(M> > *AMNI <M,N) ; 

END CASE1 ; 

ELSE 

CASE2 : IF CDISC(N) < 0 & GAMAD(M) > 0 THEN 

DO; 

MNR (N, M) = GAMA(M)*AMNR(M,N)-CN<N)*AMNI (M,N) ; 

MNI <N , M) » CN<N)*AMNR(M,N)+GAMA<M)*AMNI (M,N) ; 

END ; /■*** CASE2 **-**/ 

ELSE 

CASE3 : IF CD ISC <N) >0 & GAMAD(M) < 0 THEN 

DO ; 

MNR ( N , M > =CN ( N ) *AMNR (M , N) -GAMA ( M ) * AMN I ( M , N ) ; 

MNI <N,M)=GAMA(M)*AMNR(M,N)+CN<N)*AMNI (M.N) ; 

END ; /*** CASE 3 **■*/ 

else’ 

CASE4 : DO ; 

MNR < N , M ) *— ( CN ( N > +GAM A < M ) )*AMNI (M,N> ; 

MNI (N ,M) = (CN (N) +GAMA <M) )*AMNR(M,N) ; 

END CASE4 ; 

END LOOPMNA ; 

END LOOPAMN ; 

/** TO PRINT REAL AND IMAGINARY PART OF MATRIX MN(M,N) **/ 

/** s MNR (M , N) AND MN I ( M , N ) s M , N = -6 TO 5 **/ 

/*** ■***/ 

PUT SKIP; 

PUT SKIP ; 

PUT SKIP EDIT (REPEAT ( # , 55) ) ( X (3) , A) ; 

PUT SKIP EDIK ' PRINT VALUES OF MNREAL AND MNIMAG ') 

(X<6) , A) ; 

PUT SKIP EDIT (REPEAT < ,55) ) (X (3) ,A> ; 

MNPRINTs DO N = -5 TO 5 ; 

PUT SKIP; 

PUT SKIP ED I T ( ' N= ' , N ) <X(7),A,F(3,0>) ; 

PUT SKIP EDIT ( REPEAT ('*',10))<X(5),A>; 

PUT SKIP ; 

PUT SKIP EDIK # M' , ' MNREAL ' , 'MNIMAG' > <X<3) f A, 2 ( X (5) , A < 10) ) ) 


.01 
.02 
.03 
.04 
.05 
.06 
. 07 
.08 
.09 
. 1 

. 101 
. 102 
. 103 
. 104 
. 105 
. 106 


PUT SKIP ; 

NMPRINT : DO M = -5 TO 5; 

PUT SKIP EDIT (M, MNR (N,M) ,MNI <N,M) ) (X (2) ,F <2,0> ,2 E<15,5> > i 
END NMPRINT; 

END MNPRINT ; 

/** TO COMPUTE RIGHT HAND SIDES VECTORS FOR T(II> AND *******/ 

/** Kill) ,B2 AND B3 , REAL AND IMAGINARY PARTS -*-****/ 

/* B2 (N) = (GAMA (0) -CN (N) ) *AMN (0 , N) , N=1 ,0,-1 , -2 *****/ 

/* B3 (N ) = (GAMA < -1 ) — CN (N) ) *AMN ( - 1 , N) , N=l,0,-l,-2 *****/ 

/****•#■* *•*•***/ 


L00PB2 


: DO N = 1 TO -2 BY -1 ; 

IF CDISC(N) > O THEN 
DO ; 

B2R (N) = ( GAMA ( 0 ) -CN ( N ) ) *AMNR ( 0 , N ) ; 

B2I ( N > = (GAMA ( 0 ) — CN ( N ) )*AMNI <0,N) ; 

END ; 

ELSE 
DO ; 

B2R ( N ) =G AM A < 0 ) * AMNR ( 0 , N ) +CN < N ) *AMN I ( 0 , N ) ; 

B2I <N)=GAMA<0)*AMNI (O , N ) -CN (N) *AMNR (0 , N ) ; 

END ; END L00PB2 j 
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L00PB3: DO N=1 TO -2 BY -1; /* B3(N) FOR REAL & IMAGINARY GAMA <-!)*/ 
IF CD I SC (N) > 0 & GAMAD(-l) > 0 

THEN DO ; 

B3R<N)*(GAMA(-1>-CN<N> >*AMNR<-1 ,N) ; 

B3I ( N ) = ( GAMA ( — 1 ) — CN ( N ) )*AMNI(-1 ,N) ; 

END ; /* C(N) & GAMA (-1 > ARE BOTH REAL */ 

ELSE 

IF CD I SC (N> < 0 V. GAMAD ( -1 ) > O 
THEN DO ; 

B3R (N) —GAMA ( — 1 ) *AMNR < — 1 , N ) +CN ( N > * AMN I (-1 ,N) ; 

B3I (N) =GAMA (-1) *AMNI (-1 ,N> -CN (N) *AMNR (-1 ,N) ; 

END ; /* END OF C (N) IMAGINARY GAMA (-1 ) REAL */ 

ELSE 

IF CDISC(N) > 0 S< GAMAD (-1) < O 
THEN DO ; 

B3R (N) =- (GAMA (-1 ) *AMNI (-1 ,N> +CN (N) *AMNR (-1 ,N> > ; 

B3I (N) =GAMA<-1 > *AMNR (-1 ,N> -CN(N) *AMNI (-1 ,N> ; 

END ; /* C <N) REAL & GAMA (-1 ) IMAGINARY */ 

ELSE 
DO ; 

B3R(N)=(CN(N>-GAMA(-1) >*AMNI (-1 ,N> ? 

B3I (N) = (GAMA < — 1 ) — CN (N) >*AMNR(-1,N> ; 

END ; /* C (N) GAMA(-l) ARE BOTH IMAGINARY */ 

END L00PB3 ? /* END COMPUTING B2 AND B3 */ 

/** TO PRINT B2R(N) ,B2I (N) ,B3R(N) , AND B3I (N) , N»l f 0 f -l,-2 **/ 

/** **/ 

PUT SKIP ; 

PUT SKIP ; 

PUT SKIP EDIT (REPEAT ( ,55) ) (X(3> ,A) ; 

PUT SKIP EDIT ( ' PRINT REAL AND IMAGINARY PART OF B2 *<B3 ' ) 

(X (6) , A) ; 

PUT SKIP EDIT (REPEAT ( ' * ' , 55 ) ) ( X ( 3 ) , A ) ; 

PUT SKIP ? 

PUT SKIP ED I T ( ' N ' , ' B2REAL ' , ' B2 I MAG ' , ' B3RE AL ' , ' B3 I MAG ' > 

<X<3),A,4 (X<3) y A(10> )) ; 

PUT SKIP ; 

PRINTBs DO N s -2 TO 1 ? 

PUT SKIP EDIT (N,B2R(N> ,B2I (N) , B3R(N) , B3I (N) ) 

(X (2) , F (2,0) ,4 E ( 15 , 5) > ; 

END PRINTB; /** END PRINTING B2 AND B3 **/ 

/** PRINT AMNR(M,N) , AMNI (M,N) : M , N = -6 TO 5 **********/ 

/** **/ 

PUT SKIP? 

PUT SKIP ; 

PUT SKIP EDIT (REPEAT ( ,50) ) (X (3) ,A) ; 

PUT SKIP EDIT ( 'PRINT VALUES OF AMNREAL AND AMNIMAG ' ) 

(X (6) , A) ; 

PUT SKIP EDIT ( REPEAT < ' * ' , 55 ) ) ( X ( 3 ) , A ) ; 

PRINTAMs DO M = -5 TO 5 j 
PUT SKIP ; 

PUT SKIP EDIT ( 'M— ' ,M) (X (7) , A,F (3,0) ) ; 

PUT SKIP EDIT ( REPEAT ('*',10)) ( X ( 5 > , A ) ; 

PUT SKIP? 

PUT SKIP EDIT < ' N ' , ' AMNREAL ' , ' AMNIMAG ' ) 

( X ( 3) , A , 2 ( X (5) , A ( 10) ) ) ; 

PUT SKIP; 

PRINTAN: DO N = -5 TO 5 ? 

PUT SKIP EDIT (N,AMNR (M,N) , AMNI (M,N) ) 

< X (2) , F (2,0) ,2 E ( 15,5) ) ; 

END PRINTAN ? 

END PRINTAM ; 

/***•************************************************.*.**.*.*.*.*.*.*/ 

/* SUBROUTINE TO INPUT DATA */ 

INPUT; PROC(M); 

DCL M FIXED (4,0) ; 

GET LIST (M) ; 
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PI - 3. 14156 ; 

THETA = PI * DEGREE/ 180. ; 

SO = KO *SIN (THETA) ; 

SM(M) « SO + M*2*PI/D ; 

K02 = (D2/D) *K0**2* (E1-E2) ; 

U1(M) = SQRT(SM(M) **2+K02>. ; 

IF FLAG= 1 THEN DO; 

PUT SKIP (2) ; 

PUT SKIP<2> ; 

PUT PAGE EDIT ( ' NEWTON ' ' S METHOD ' > < X ( 10) , A) ; 

PUT SKIP (2) ; 

PUT SKIP (2) EDIT('***TE MODE ***')< X < 1 0 >, A ) ; 

PUT SKIP (2) ; 

PUT SKIP <3> EDIT ( ' INITIAL VALUES ' ) ( X ( 10 ) , A) ; 

PUT SKIP (2); 

PUT SK IP ( 2 ) EDIT ( ' M= ' , M ) ( X ( 5 > , A , F ( 4 , 0 ) ) ; 

PUT SKIP EDIT ( 'U1 <M) « ' ,U1 (M) , 'EPSILQN= EPSILON, 

' KM AX— ' , KMAX ) (X (5) ,2 (A ,E ( 12 , 5) , X (5) ) , A , F (3 , 0) » 


339 . 


PUT SKIP (5 ) EDIT< 'U1 ' , 'U2' , 'FUJI) ' , 'F' ' <U1> ' , 

340. 


' MJ2-U1 ! ' , * COUNT ' ) (X (3) ,A,X ( 12) , A, X ( 10) , A, 

341 • 


X (7) , A , X (8 ) , A, X (2) , A) ; 

342. 


END; 

343. 


RETURN; 

344. 


END INPUT; 

345. 


/* */ 

346. 


/* SUBROUTINE TO PERFORM CALCULATION */ 

347. 


CALCULATE: PROC(M); 

348. 


DCL M FIXED (4,0) ; 

349. 


U2(M)=U1 (M)-F(Ul (M) , M) /FPR I ME ( U 1 (M) ,M) ; 

350 • 


K=K+ 1 ; 

351. 


RETURN; 

352. 


END CALCULATE; 

353. 


/* */ 

354. 


/* SUBROUTINE TO PRINT TABLE */ 

355. 


PRINT: PROC(M) ; 

356. 


DCL M FIXED (4 ,0> ; 

357. 


PUT SKIP EDIT ( U 1 ( M ) , U2 (M) ,F(U1 (M) ,M) , FPR I ME ( U 1 ( M ) , M ) * 

358. 


ABS(U2(M)-U1 (M) ) ,K> 

359. 


(5 (E ( 12, 5) ,X(1> ) , F (2,0) ) ; 

360 • 


RETURN; 

361. 


END PRINT; 

362. 


/* */ 

363 • 


/* SUBROUTINE TO PRINT FINAL RESULTS */ 

364. 


OUTPUT: PROC(M); 

365. 


DCL M FIXED (4,0) ; 

366. 


PUT SKIP (5) EDIT ( 'APPROXIMATE ROOT U2= ' ,U2(M), 

367. 


' F <U2) « ' ,F(U2(M) ,M) ) (A,E(14,7) ,X(5) ,A,E(14,7) ) ; 

368. 


RETURN; 

369. 


END OUTPUT; 

369.01 

/* 

TO COMPUTE SCATTRING MATRIX SCAT */ 

369.02 

/* 

CONSTRUCT COEFFICIENT MATRIX AA FROM MATRIX MN(N,M) */ 

369.03 


MM= 8; 

369.04 


NN= 8 ; 

369.041 


AA 11: DO I = 1 TO NN/2 ; 

369.042 


AAJ1: DO J = 1 TO MM/2; 

369.043 


A A ( I , J ) =MNR ( 2- I , 2-J > ; 

369.044 


AA ( I , J+MM/2) =-MNI (2-1 ,2-J) ; 

369.045 


END AA J 1 ; 

369.046 


END AA I 1 

369 . 05 


P I VOT— 1.0 ; 

369 . 06 


RS= 3 ; 

369. 11 


AA 12: DO I = NN/2+1 TO NN ; 

369. 12 


AAJ2: DO J = 1 TO MM/2 ; 

369. 13 


AA ( I . J ) =MNI < NN/2+2— I .2-J) s 


320. 

321. 

“TOO 

324. 

325. 

326. 

327. 

328. 

329. 

330 . 

331 . 


334 . 

335. 

336. 

337. 
338 • 
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369. 14 


AA(I , J+MM/2) ~MNR <NN/2+2- 1 ,2-J) ; 

369. 15 


END AAJ2 ; 

369 . 1 6 


END AAI2 ; 

369.21 


RHS1 : DO I = 1 ,3 TO NN ; 

369.22 


AA ( I , MM+ 1 ) =0 ; 

369.23 


END RHS1 ; 

369.24 


AA<2,MM+1)=2*CN(0> ; 

369.25 

/* 

END OF RHS1 COLUMN ' */ 

369.29 


RHS23R : DO I » 1 TO NN/2 ; 

369.31 


AA(I ,MM+2)=B2R(2-I) ; 

369.32 


AA(I ,MM+3> =B3R <2-1 > ; 

369.33 


END RHS23R ; 

369. 34 

/* END OF RHS2 COLUMN OF AUGMENTED MATRIX IN GAUSS PROC. */ 

369.35 


RHS23I: DO I = NN/2+1 TO NN ; 

369 . 36 


AA ( I , MM+2 ) =B2 1 ( MM/2+2- 1 ) ; 

369.37 


AA ( I , MM+3) =B3 I (MM/2+2-1 > ; 

369.372 


END RHS23I ; 

369.38 

/* 

END OF RHS3 COLUMN OF AUGMENTED MATRIX IN GAUSS PROC. •*/ 

369.39 


CALL GAUSS(AA, MM, NN, PIVOT, RS,X) ; 

369.4 

/* 

INVOKE GAUSS ELIMINATION TO COMPUTE SCATTERING MATRIX */ 

400. 

GAUSS s 

PROC ( AA , M , N , P I VOT , RS , X ) ; 

404. 


/* */ 

407. 


/* GAUSSIAN ELIMINATION WITH OR WITHOUT PIVOTING. *-/ 

408. 


/* ANSWERS ARE THEN SUBSTITUTED BACK INTO THE */ 

409. 


/* ORIGINAL EQS. WITH MULTIPLE RHS VECTORS. */ 

410. 


/* */ 

456. 


DCL (AA (*,*) , X (*,*) ) FLOAT 5 

457. 


DCL(M,N, PIVOT, RS) FIXED<5,0> ; 

459. 


STARTs BEGIN; 

460. 


DCL (AC (M , N+RS) ,BB(M,N+RS) ,HOLD(N+RS) , SUM) FLOAT (6 ) 

461. 


' XX (RS , N> FLOAT (6) INIT ( <RS*N) 0) ; 

462. 


/* INPUT AUGMENTED MATRIX */ 

463 . 


CALL INPUT1; 

464. 


/* CONVERT TO UPPER TRIANGULAR MATRIX */ 

465. 


CALL UPTRI ; 

466. 


/* BACK SUBSTITUTE */ 

467. 


CALL BACKSUB; 

46Q. 


CALL OUTPUT 1; 

469. 


/* PUT ANSWERS BACK IN ORIGINAL EQUATIONS */ 

470. 


CALL TEST!; 

477. 


/* */ 

478. 


/* SUBROUTINE TO INPUT AUGMENTED MATRIX */ 

479. 


INPUT Is PROC; 

480. 


PUT PAGE EDIT < 'GAUSSIAN ELIM INAT I ON ')( X ( 28) , A) ; 

481. 


IF PIV0T=1 THEN PUT SKIP EDIT('WITH PIVOTING') 

482. 


(X <31> , A) ; 

483. 


ELSE PUT SKIP EDIT( 'WITHOUT PI VOTING ')( X (30) , A> ; 

484. 


PUT SKIP ED I T ( ' FOR ' f M,' BY ',N,' MATRIX') 

485. 


( X (29) ,A,F(2,0) ,A,F(2,0> ,A) ; 

486. 


PUT SKIP EDIT < 'WITH' ,RS, 'RIGHT HAND SIDES') 

487. 


(X (29) , A , F < 3 , 0 > ,X (2> ,A) ; 

488. 


PUT SK I P ( 5 ) ; 

489. 


DO 1=1 TO M; 

490. 


DO J=1 TO N+RS; 

491. 


AC < I , J > =AA ( I , J > ; 

492. 


PUT EDIT (AA(I , J) ) (X < 1) ,F(8,3) ) ; 

493. 


BB ( I , J ) =AC ( I , J > ; 

494. 


END; 

495. 


PUT SKIP; 

496. 


END; 

497. 


RETURN: 

498. 


END INPUT1; 

499 . 


/* */ 

500 . 


/* SUBROUTINE TO PRINT MATRIX */ 

501. 


PRINT: PROC; 

502. 


PUT SKIP (5): 


52 


503. 
504 « 

505 . 

506. 
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510. 
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515. 
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520. 
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529. 
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541. 
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54o . 

544. 

545. 
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548. 
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552. 

553. 

554. 

555. 

556. 

557. 

558. 

559. 

560. 

561. 
561. 1 

562. 

563. 

564. 

565. 

566. 

567. 


DO 1=1 TO M; 

DO J®1 TO N+RS; 

PUT ED I T ( AC ( I , J ) ) <X(i) ,F<8,3>); 

END; 

PUT SKIP; 

END; 

RETURN; 

END PRINT; 

/* */ 
/* SUBROUT. CONVERTS MATRIX TO UPPER TRIANGULAR */ 
UPTRIs PROC; 

DO K=1 TQ M-l ; 

IF PI V0T=1 THEN CALL PIV0T1; 

DO I«K+1 TO M; 

RATIO = AC ( I , K ) /AC ( K , K ) ? 

DO TO N+RS; 

AC < I , J > =AC ( I , J ) —RAT 10* AC <K , J > ; 

END; 

END; 

END; 

RETURN; 

END UPTRI ; 

/* SUBROUTINE TO USE PIVOTING */ 

PIVOTls PROC; 

P=K ; 

DO I-K+l TO M; 

IF ABS (AC (P , K) ) < ABS ( AC < I , K ) > THEN P = I; 

END; 

IF P' Vs =K THEN DO; 

DO J=1 TO N+RS; 

HOLD ( J ) =AC ( K , J ) ; 

AC < K , J ) =AC < P , J ) ; 

AC ( P , J ) =HOLD < J ) ; 

END; 

END; 

RETURN; 

END PIVOT 1 ; 

/* - */ 
/* SUBROUTINE TO BACK SUBSTITUTE */ 

BACKSUB: PROC; 

DO K * 1 TO RS ; 

DO I=N TO 1 BY(-l) ; 

SUM=0; 

DO J=I TO M; 

SUM*SUM+X X ( K , J ) *AC < I , J ) ; 

END; 

X X <K , I ) = < AC ( I , N+K) -SUM ) /AC (1,1) ; 

END; 

END ; 

RETURN; 

END BACKSUB; 

/* */ 
/* SUBROUTINE TO PRINT ANSWERS */ 

OUTPUT 1: PROC; 

PUT SKIP (5) EDIT( 'ANSWERS' ) (X (34) , A) ; 

DO J « 1 TO.RS ; 

PUT SKIP EDIT ( ' SET ' , J) < X < 20 > , A , F ( 3 , 0 > ) ; 

PUT SKIP; 

DO 1*1 TO N; 

X (J, I)=XX (J , I ) ; 

PUT SKIP EDIT < ' X < ' , J , * T ' , I , ' > = , X X ( J , I ) ) 

(X (28) , A , F ( 2 , 0 ) , A , F < 2 , 0 ) , A , F ( 9 , 6 ) ) ; 

END; 

END ; 

RETURN; 

END OUTPUT 1: 
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/* */' 
/* SUBROUTINE TO PUT ANSWERS BACK IN ORIGINAL */ 
/* EQUATIONS */ 

TEST 1 : PROC; 

PUT SKIP (5) EDIT 

('ANSWERS PUT IN ORIGINAL EQUATIONS 'M X <21 >, A) 
DO K =* 1 TO RS ; 

PUT SKIP EDIT ('SET ' „ K) ( X ( 29 ) , A , P ( 2 , 0 ) ) ; 

PUT SKIP? 

DO 1=1 TO M; 

PUT SKIP; 

SUM=0 ; 

PUT SKIP EDIT< ' '>(X(1>,A>; 

DO J=1 TO N; 

SUM=SUM+BB ( I , J ) *XX (K , J ) ; 

PUT EDIT <BB ( I , J > , ' X ( ' , J , ') ') 

(P (9,3) , A , F ( 1 , 0 ) , A > ; 

IF J < N THEN PUT EDIT('+ ' ) (A) ; 

ELSE IF J=N THEN PUT EDIT('= ' > (A) ; 

END; 

PUT EDIT (SUM) (F (9 ,3) ) ; 

END; 

END ; 

RETURN; 

END TEST1 ; 

END START ; 

END GAUSS ; 

/* */ 
CALL LOOPT; /* PROCEDURE TO COMPUTE T VECTORS */ 

/* TO PRINT OUT T VECTORS PROCEDRE CALL */ 

CALL PR I NTT; 

CALL INNER ; /* PROCEDRUE TO COMPUTER DOT PRODUCT */ 

CALL SCATTER; /* PROCEDURE TO COMPUTE SCATTERING MATRIX*/ 
CALL RTL ; /* PROCEDURE TO COMPUTE RTL MATRIX */ 

/* PROCEDURE TO COMPUTE T VECTORS */ 

/* TR AND TI ARE REAL AND IMAGINARY PART OF T VECTORS */ 

/* */ 
LOOPT: PROC; 

DO K 33 1 TO 3 ; * 

DO N = 1 TO -2 BY -1 ; 

TR(K,N)*X<K,2-N> ; 

TI <K,N)=X (K,6-N) ; 

END; 

END; 

RETURN; 

END LOOPT ; 

/* PRINT T VECTORS, TREAL AND TIMAGINARY */ 

/* */ 

PR I NTT; PROC; 

PUT SKIP; 

PUT SKIP EDIT (REPEAT ( ' ,55) ) (X (3) , A) ; 

PUT SKIP EDIT ( ' PRINT T VECTORS FOR N - 1 ,0,-1 ,-2 ') 

(X (7) ,A> ; 

PUT SKIP ED I T ( REPEAT ( ' * ' , 55 ) ) ( X < 3 > , A ) ; 

PUT SKIP ; 

DO K = 1 TO 3 ; 

PUT SKIP; 

PUT SKIP EDIT ( ' K> ' , K ) ( X ( 1 0 ) , A , F ( 3 , 0 > > ; 

PUT SKIP EDIT (REPEAT ( '*' , 10) ) (X (5) ,A> ; 

PUT SKIP; 

PUT SKIP EDIT ( 'N' , 'TREAL' , 'TIMAG' ) (X (7) ,A,2 (X (4) ,A(10) ) ) 
PUT SKIP; 

DO N = -2 TO 1 ; 

PUT SKIP EDIT ( N , TR ( K , N ) , T I ( K , N ) ) ( X ( 5 ) , F ( 3 , O > ,2 E(15,5)>; 
PUT SKIP; 

END: 
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END; 


END PRINTT; 

/* PROCEDURE TO COMPUTE INNER PRODUCT OF AMN AND T MATRICES*/ 
/* •*/ 


INNER : PROC; 

DO K = 1 TO 3 ; 

ROR(K) = 0; 

RQI (K) = O' ; 

DO N = -2 TO 1 ; 

ROR ( K ) =ROR ( K ) +AMNR ( N , O ) *TR ( K , N ) -AMN I (N , 0) *TI ( K , N) ; 
ROI ( K) =ROI (K> +AMNI < N , 0 ) *TR ( K , N > +AMNR ( N , 0 ) *T I CK.N) ; 
END ; 

END; 

ROR(1)=ROR<1)-1; 

ROR ( 2 ) =ROR ( 2 ) + AMNR (0,0) ; 

ROI (2) =ROI (2)+AMNI (0,0) ; 

ROR ( 3 ) =ROR ( 3 ) +AMNR ( — 1 ,0) ; 

ROI (3) =ROI (3) +AMNI (-1,0) ; 

/* PRINT INNER PRODUCT OF AMN ( N , 0 ) AND T ( K , N ) */ 

/* */ 


PUT SKIP; 

PUT SKIP EDIT (REPEAT ( ,55) ) (X(7) ,A) ; 

PUT SKIP EDIT ( 'RO, INNER PRODUCT OF A & T')(X(5),A>; 

PUT SKIP EDIT (REPEAT ( '*',55)) (X(7) ,A) ; 

PRINTRO: DO I = 1 TO 3 ; 

PUT SKIP EDIT ( ' ROREAL ( ' , I , ' ) = ' , ROR ( I ) , ' ROIMAG ( ' , I , ' ) = ' , 
ROI (I) ) (X(10) ,2 ( X (2) , A , F (2,0) ,A,E(12,5))); 

END PRINTRO ; 

END INNER ; 

/* END ON INNER PRODUCT PROCEDURE */ 

/* PROCEDURE TO COMPUTE SCATTERING MATRIX */ 

SCATTER: PROC; 

L00PSC1 s DO j = -1 TO 1 ; 

SCATR ( 1 , J ) =ROR ( 2- J ) ; 

SCAT I ( 1 , J ) =R0 I <2-J ) ; 

END L00PSC1 ; 

LQGPSC2: DO I * 0,-1 ; 

DO J = -1 TO 1 ; 

SCATR < I , J ) =TR ( 2- J , I > ; 

SCATI ( I , J) -TI (2— J , I ) ; 

END ; 

END L00PSC2; 

/* PRINT SCATTERING MATRX SCATREAL , SCATIMAG PARTS */ 


/* 


*/ 


PUT SKIP; 

PUT SKIP EDIT (REPEAT ( ' * ' , 55) ) ( * (3) , A) ; 

PUT SKIP EDIT ('PRINT SCATTERING MATRIX SCAT') 

(X (10) , A) ; 

PUT SKIP EDIT (REPEAT ( ,55) ) (X<3) ,A> ; 

PUT SKIP; 

DO K * -1 TO 1 ; 

PUT SKIP; 

PUT SKIP EDIT ( ' K= ' , K ) ( X ( 10) , A , F (2 , 0) ) ; 

PUT SKIP EDIT (REPEAT ( ,10) ) (X(5) ,A) ; 

PUT SKIP; 

DO N — -i TO 1 ; 

PUT SKIP EDIT ( ' SREAL ( ' ,N, ' , ' ,K, ' ) = ' , SCATR (N,K) , 

' SI MAG ( ' ,N, ' , ' , K , ' ) = ' , SCATI ( N , K ) ) 
(X(5) ,2 (X(2) , A , F ( 2 , O ) , A , F ( 2 , 0 ) , A , E ( 1 2 , 5 ) ) ); 
END ; 

END; 

RETURN; 

END SCATTER; 

/* COMPUTE COEFFICEINT MATRIX RTL, AND THE RHS VECTOR**/ 
/* OF THE SYSTEM OF EQUATIONS TO SOLVE R,R',T(0), AND */ 
/* T ( — 1 ) */ 
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RTL : PROC; 

RTLR (1 , 1 ) =1 . 0 ; 

RTLR (2 , 2 ) =1 . 0 ; 

DO I = 2 TO 4 ; 

RTLR ( I , 1 ) =0 ; 

END ; ' /* RTLR (2,1) =RTLR (3,1) =RTLR <4.1) =0 */ 

DO I = 1, 3 TO 4 ; 

RTLR (1,2) “0 3 /* RTLR (1,2) =RTLR (3,2) =RTLR (4,2) =0 */ 

END ; 

DO I = 1 TO 4 ; 

RTL I ( 1 , 1 ) =0 ; ' 

RTL I < 1 , 2 ) =0 ; 

END ; /* RTL I ( T, 1 ) =RTLI (1,2) =0 , 1 = 1 ,2,3,4 •*/ 

/* END OF DEFINING RTL MATRIX INITI ALI ZATION */ 

GAMA0H=GAMA (0) *H ; 

GAMA_ 1 H=GAMA ( - 1 ) *H ; 

GAMA_ 1 2H=2*GAMA_ 1 H ; 

GAMA02H=2*GAMA0H ; 

GAMAO_ 1 H= < GAMA < 0 ) +GAMA ( - 1 ) ) *H ; 

/* TO COMPUTE RTL (1,3) */ 

/* REAL & IMAGINARY PARTS OF S ( 1 , 0) *S (0 , 0) & S < 1 , -1 ) *S (-1 , 0 > */ 
DO I = O TO -1 BY -1 ; 

AB ( 1 - 1 ) =SCATR (1,1) *SCATR (1,0) -SCAT 1(1,1) *SCAT 1(1,0) ; 

CD ( 1-1 ) -SCATR (1,1) *SCAT I (1,0) +SCATI (1,1) *SCATR (1,0) ; 

END ; /* END OF DEFINING AB ( 1 ) , AB (2) , CD ( 1 ) , & CD(2) */ 

/* RTL (1,3) */ 

if gamad(-i) > o Then 

DO ; /* RTL (1,3) , WHEN GAMA (-1 ) IS REAL */ 

RTLR (1,3)=- (AB ( 1 ) *COS ( GAMA02H ) -CD ( 1 ) *S I N ( GAMA02H ) 

+ AB (2) *COS (GAMA0_1H) —CD (2) *SIN (GAMA0_1H) ) ; 

RTL I ( 1 , 3 ) =- ( AB ( 1 ) *S I N ( GAMA02H ) +CD ( 1 ) *COS ( GAM A02H > 

+ AB(2)*SIN(GAMAO_1H)+CD(2)*COS<GAMAO_1H) ) ; 

END ; 

ELSE 
DO ; 

RTLR (1,3) =CD < 1 ) *S I N ( GAMA02H ) - AB ( 1 ) *COS ( GAMA02H ) 

+EXR ( -GAMA (- 1 ) ) * (CD (2) *SIN (GAMAOH) -AB (2) * 

COS (GAMAOH) ) ; 

RTL I ( 1 , 3 ) =- ( AB ( 1 ) *S I N ( GAMA02H ) +CD ( 1 ) *COS ( GAMA02H ) 

+EXP (—GAMA (—1 ) )* (CD (2) *SIN (GAMAOH) +AB (2) *C05 (GAMAOH) ) ) ; 

END; 

/* TO COMPUTE RTL (1,4) */ 

/* REAL & IMAGINARY PARTS OF S < 1 f 0) *S <0,-1 ) & S < 1 , -1 ) *S (-1 , - 1 ) */ 
DO I = O TO -1 BY -1 ; 

AB (3-1 ) =SCATR (1,1) *SCATR (1,-1) -SCATI (1,1) *SCATI (1,-1); 

CD ( 3— I ) =SCATR (1,1) *SCAT I (1,-1) +SCATI ( 1 , I ) *SCATR (1,-1); 

END ; 

/* RTL (1,4) */ 

IF GAMAD(-l) > 0 THEN 

DO ; /* RTL (1,4) WHEN GAMA(-l) IS REAL */ 

RTLR < 1 , 4 ) =- ( AB ( 3 ) *COS ( GAMAO_ 1 H ) -CD (3) *SIN ( GAMAO_ 1 H ) 

+AB<4>*C0S(GAMA_12H)H:D(4)*SIN<GAMA_12H) ) ; 

RTL I ( 1 , 4 ) =- ( AB ( 3 ) *S I N ( GAMAO__ 1 H > +CD ( 3 ) *COS (GAMA 0 _ 1 H) 

+AB (4) *SIN (GAMA I2H) +CD (4) *COS (GAMA_ 12H ) ) ; 

END ; “ ’ 

ELSE 

DO ; /* RTL (1,4) WHEN GAMA (-1 ) IS IMAGINARY */ 

RTLR ( 1 , 4) =EXP (-GAMA (- 1 ) ) * (CD (3) *S I N (GAMAOH ) -AB (3) 

*COS( GAMAOH) ) -AB ( 4 ) *EXP ( ~GANA_ 1 2H ) ; 

RTL 1(1,4) “--EXP ( -GAMA ( - 1 ) ) * ( AB ( 3 > *S I N < GAMAOH ) +CD ( 3 ) 

★CDS ( GAMAOH > ) -CD < 4 > *EXP < -GAMA_ 1 2H ) ; 

END ; 

/* TO COMPUTE RTL (2,3) AND RTL(2,4) **/ 

RTLR ( 2 , 3 ) =- ( SCATR (1,0) *COS ( GAMAOH ) -SCATI ( 1 , 0 ) *S IN (GAMAOH) ) ; 

RTL I (2.3) =-( SCATR (1.0) *S I N (GAMAOH) +SC ATI ( 1 . 0 ) *COS ( GAMAOH ) ) .- 
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756. IF GAMAD(-l) > 0 THEN 

757. DO ; /* RTL<2,4> WHEN GAMA(-l) IS REAL */ 

758. RTLR (2, 4) — - ( SC ATR (1,-1) *CQS ( GAM A_ 1 H ) -SCAT I ( 1 ,-l ) * 

759. SIN (GAMA_1H) ) ; 

760. RTL I (2,4) =- (SCATR (1,-1) *SIN (GAMA_1H) +SCATI < 1 , - 1 ) *COS < GAMA_ 1 H ) ) 
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END ; 

ELSE 

DO ; /* RTL (2, 4) WHEN GAMA (-1 ) IS IMAGINARY */ 

RTLR (2,4) =-EXP ( -GAMA_ 1 H > *SCATR (1,-1); 

RTLI (2,4) --EXP (-GAMA 1H)*SCATI (1,-1) ; 

END ; 

/* TO COMPUTE RTL (3,3) */ 

/* REAL & IMAGINARY PARTS OF S(0,0)**2 & S <0,-1 > *S (-1 ,0) */ 

DO I = 0 TO -1 BY -1 ; 

AB ( 5— I ) —SCATR < 0 , I ) *SCATR (1,0) —SCAT I (0, I)*SCATI (1,0) ; 

CD ( 5— I ) =SCATR (0,1) *SCAT I (1,0) +SCATI (0, I ) *SCATR (1,0); 

END ; /* END OF COMPUTING REAL & IMAGINARY PARTS OF SCAT */ 

/* RTL (3,3) */ 

IF GAMAD(-l) > O THEN 

DO ; /* RTL (3, 3) WHEN GAMA (-1 ) IS REAL */ 

RTLR ( 3 , 3 ) =- AB ( 5 ) *COS ( G AM A02H ) +CD ( 5 ) *S I N ( GAMA02H ) 

-AB (6) *COS ( GAM A0_ 1 H ) +CD (6) *SIN (GAMA0_1H> +1 ; 

RTL I ( 3 , 3 ) — ( AB < 5 ) *S I N ( G AM A02H ) +CD ( 5 ) *COS ( G AM A02H ) 

+ AB ( 6 ) *S I N ( G AM A0_ 1 H ) +CD ( 6 ) *COS ( GAMAO_ 1 H ) ) ; 

END ; 

ELSE 

DO ; /* RTL (3,3) WHEN GAMA(-l) IS IMAGINARY */ 

RTLR (3,3) -1-AB (5) *COS (GAMA02H) +CD (5) *SIN (GAMA02H) -EXP (-GAMA_1H 

* ( AB (6) *COS (GAMAOH) -CD (6) *SIN (GAMAOH) ) ; 

RTL 1(3,3) *-AB (5) *S IN ( GAMA02H ) -CD ( 5 ) *COS ( GAMA02H ) 

-EX P ( ~GAMA_ 1 H > * ( AB ( 6 ) *S I N ( GAMAOH ) +CD (6) *COS < GAMAOH ) > ; 

END ; 

/* TO COMPUTE RTL (3,4) */ 

/* REAL & IMAGINARY OARTS OF S (0 , 0) *S (0 , -1 ) & S (0 , -1 ) *S (-1 , -1 ) */ 
DO I = 0 TO -1 BY -1 ; 

AB ( 7- I ) =SCATR (0,1) *SCATR ( I , — 1 ) -SCATI (0,I)*SCATI (1,-1) ; 

CD ( 7— I ) =SCATR ( 0 , I)*SCATI (I,-1)+SCATI (0,1) *SCATR (1,-1) ; 

END ; 

/* RTL (3,4) */ 

IF GAMAD(-l) > 0 THEN 

DO /* RTL (3 , 4) WHEN GAMA (-1 ) IS REAL */ 

RTLR ( 3 , 4 > *- < AB ( 7 ) *COS ( GAMAO_ 1 H ) -CD (7) ♦SIN ( GAMAO_ 1 H ) 

+AB(8)*C0S(GAMA_12H)-CD(B)*SIN(GAMA_12H> ) ; 

RTLI (3,4) «- (AB (7) *SIN (GAMAO^IH) +CD (7) *COS ( GAMAO_ 1 H ) 

+AB(8)*SIN(GAMA_.12H)+CD(8)*C0S(GAMA_12H> ) ; 

END ; 

ELSE 

DO ; /* RTL (3,4) WHEN GAMA(-l) IS IMAGINARY */ 

RTLR (3,4) =*EXP ( -GAMA_ 1 H ) * ( CD ( 7 ) *S I N ( GAMAOH ) -AB ( 7 ) *COS ( GAMAOH ) ) 
-EXP ( -GAMA_ 1 2H ) *AB ( 8 ) ; 

RTLI (3,4) a-EXP (-GAMA_1H) * (AB (7) *SIN (GAMAOH) +CD (7) * 

COS ( GAMAOH ) ) -CD ( 8 ) *EXP ( ~GAMA_ 1 2H ) ; 

END ; 

/* TO COMPUTE RTL (4,3) */ 

/* REAL & IMAGINARY PARTS OF S ( - 1 , 0 ) *S ( 0 , 0 > & S ( - 1 , - 1 > *S ( - 1 , 0 ) */ 
DO I = 0 TO -1 BY -1 ; 

AB(9-I) =SCATR(-1 , I ) *SCATR (1,0) -SCAT I (-1,1) *SCATI ( I ,0) ; 

CD (9-1) -SCATR (-1 , I ) *SCATI ( 1,0) +SCATI (-1,1) *SCATR < I ,0) ; 

END ; 

/* RTL (4,3) #/ 

IF GAMAD(-l) > O THEN 

DO ; /* RTL(4,3) WHEN GAMA (-1 ) IS REAL */ 

RTLR (4,3) =»- (AB (9) *COS ( GAMA02H ) -CD (9) *SIN ( GAMA02H ) 

+AB (10) *COS (GAMAO 1H) -CD (10) *-SIN (GAM AO 1H) ) s 
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RTL I (4,3) =- (AB (9) *SIN < GAMA02H) +CD (9) *C0S (GAMA02H) 

+AB (10) *S I N ( GAMAO 1 H ) +CD (10) *COS ( GAMAO 1 H ) ) ; 

END ; 

ELSE 

DO ; /* RTL (4,3) WHEN GAMA (-1 ) IS IMAGINARY */ 

RTLR (4,3) =— AB ( 9 ) *COS ( GAMA02H ) +CD < 9 > *S I N ( GAMA02H ) 

-EXP ( -GAMA_ 1 H ) * ( AB ( 1 0 ) *COS ( GAMAOH ) -CD < 1 0 ) *S I N < GAMAOH > ) ; 
RTL 1(4,3) =-AB ( 9 > *S I N (GAMA02H ) -CD ( 9 ) *COS ( GAMA02H > 

-EXP (-GAMA_ 1 H) * ( AB < 10) *SIN (GAMAOH) +CD ( 10) *COS (GAMAOH) ) ; 
END ; 

/* TO COMPUTE RTL (4,4) */ 

/* REAL & IMAGINARY PARTS OF S (-1 ,0) *S (0,-1 ) Z< S(-l,-l>**2 */ 
DO I « 0 TO -1 BY -1 ; 

AB(ll-I) =SCATR(-1 , I) *SCATR( I ,-l>-SCATI (-1 , I) *SCATI (1,-1); 
CD (11 — I ) =SCATR (“1 , I)*SCATI (1,-1) +SCATI (-1 , I ) *SCATR (1,-1); 
END ; 

/* RTL (4,4) */ 

IF GAMAD(-l) > 0 THEN 

DO 5 /* RTL (4,4) WHEN GAMA (-1 ) IS REAL */ 

RTLR (4 , 4) =- (AB ( 1 1 ) *COS (GAMA0_1H) -CD (11) *S IN (GAMAO_ 1 H) 

+AB ( 12) *CQS (GAMA_12H> -CD ( 12) *SIN (GAMA_12H) ) +1 ; 
RTLI (4,4)=- (AB (11) *S I N ( GAMAO_ 1 H ) +CD (11) *COS (GAMA0_1H) 

+AB ( 12) *S IN (GAMA_12H> +CD ( 12) *COS (GAMA_i 2H) ) ; 

END ; 

ELSE 

DO ; /* RTL (4,4) WHEN GAMA ( - 1 ) IS IMAGINARY */ 

RTLR (4 , 4) =1-EXP (-GAMA_ 1H) * ( AB ( 1 1 ) *COS (GAMAOH) -CD ( 1 1 ) * 

SIN (GAMAOH) ) —EXP ( -GAMA_ 1 2H ) *AB (12) ; 

RTLI (4,4) =— EXP (-GAMA^IHT * (AB ( 1 1 ) *SIN (GAMAOH) +CD ( 1 1 ) * 

COS (GAMAOH) ) -EXF*7-GAMA_12H) *CD (12) ; 

END ; 

/* END OF COMPUTING COEFFICIENT MATRIX RTL FOR SOLVING */ 

/* R AND RPR I ME, T(0) AND T(-l) IN A SYSTEM OF EQUATIONS */ 

/* TO PRINT ELEMENTS OF RTL MATRIX */ 

PUT SKIP ; 

PUT SKIP EDIT (REPEAT ( ' * ' , 55) ) (X (3) , A) ; 

PUT SKIP EDIT( 'PRINT COEFFICIENT MATRIX RTL ' ) ( X ( 10 ) , A ) ; 

PUT SKIP EDIT (REPEAT ( ,55) ) <X(3> ,A) ; 

PUT SKIP ; 

DO I s 1 TO 4 ; 

PUT SKIP; 

PUT SKIP ED I T ( ' I = ' , I ) ( X ( 1 0 > , A , F ( 2 , 0 ) ) ; 

PUT SKIP EDIT (REPEAT ( , 10) ) (X (5) ,A) ; 

PUT SKIP ; 

DO J » 1 TO 4 ; 

PUT SKIP EDIT ( ' RTLR (',I,',',J,')= # , RTLR ( I , J ) , ' RTLI ( ' , I , ' , ' , 
J, ')=' , RTLI (I,J) ) 

( X (5) ,2 (X(2) , A , F ( 2 , 0 ) ,A,F(2,0) ,A,E(12,5) ) > ; 
END ; /* END OF PRINTING RTL(I,J) FOR I,J = 1 TO 4 */ 

END ; /* END OF PRINTING RTL(I,J) FOR I,J = 1 TO 4 */ 

END RTL ; 

/* TO COMPUTE R , RPR I ME , T ( 0 ) , AND T(-l) */ 

MM= 8 ; 

NN= 8 ; 

PI VOT= 1 ; 

RS= 1 ; 

AAII1: DO 1=1 TO NN/2; 

AAJJls DO J= 1 TO MM/2; 

AA ( I , J ) =RTLR ( I , J) ; 

AA(I , J+NN/2)=-RTLI ( I , J > ; 

END AAJJ1 ; 

END AA III; 

AAII2: DO 1= NN/2+1 TO NN ; 

AAJJ2: DO J= 1 TO MM/2 $ 

AA ( I , J ) =RTL I (I -NN/2, J) ; 

A A ( I . J+NN/2) =RTLR ( I -NN/2 . J ) : 


58 



*/ 

*/ 


918. 

919. 

920. 

921. 

922. 

923. 

924. 

925. 

926. 

927. 

928. 

929. 

930. 
1042. 

1044. 

1045. 

1046. 

1047. 

1048. 

1049. 

1050. 

1051. 

1052. 

1053. 

1054. 

1055. 

1056. 

1057. 

1058. 

1059. 

1060. 
1061. 
1062. 

1063. 

1064. 

1065. 

1066. 

1067. 

1068. 

1069. 

1070. 

1071. 

1072. 

1073. 

1074. 

1075. 

1076. 

1077. 

1078. 

1079. 

1080. 
1110 . 
1111 . 
1112 . 

1113. 

1114. 

1115. 

1116. 


END AAJJ2; 

END AAII2; 

/* RIHGT HAND VECTOR FOR THE SYSTEM OF EQUATIONS 
/* TO SOLVE FOR R, RPR I ME, T(0),«c T<-1) 

AA (1 , MM+ 1 ) —SCATR (1,1); 

AA ( 5 , MM+ 1 ) —SCAT I (1,1); 

AA (2, MM+1 ) =0 ; 

A A (6, MM+1 ) “0 ; 

RHSCAT : DO I = 3 TO NN/2; 

AA ( I , MM+ 1 ) —SCATR ( 3— I , 1) : 

AA ( I+NN/2 , MM+ 1 ) -SCAT I (3-1 ,1) ; 

END RHSCAT ; 

CALL GAUSS (AA, MM, NN, PIVOT, RS, X) ; 

/* */ 
RREAL—X (1,1); /* REAL PART QF R */ 

RIMAG— X (1,5) ; /* IMAGINARY PART OF R */ 

RPRIMER=X (1,2); /* REAL PART OF R PRIME */ 

RPR I ME I = X ( 1 ,6) ; /* IMAGINARY PART OF R PRIME */ 

TOR » X(l,3>; /* REAL PART OF (1,0) */ 

TOI = X < 1 ,7) ; /* IMAGINARY PART OF (1,0) */ 

TMINUS1R=X (1,4); /* REAL PART OF (1,-1) */ 

TMINUS1 1-X ( 1 , 8) ; /* IMAGINARY PART OF (1,-1) */ 

RABS=SQRT(RREAL**2+RIMAG**2) ; 

RPHASE-ATAND (RIMAG/RREAL) ; /* PHASE ANGLE OF R */ 
RF’RABS-SQRT ( RPR I MER**2+RPR I ME I **2 ) ; 

RPRPHASE^AT AND ( RPR I MER/RPR I ME I ) ; /*PHASE OF RPR I ME */ 

/* TO PRINT R */ 

PUT SKIP; PUT SKIP ; 

PUT SKIP EDIT ('REAL OF R =' ,RREAL IMAGINARY OF R= ' , RIMAG) 
(X (3) ,2 (X (2) , A,E ( 12 ,5) ) ) ; 

PUT SKIP ; 

PUT SKIP ED I T ( ' ABS ( R ) = ' , RABS , ' PHASE ( R ) = ' , RPH ASE ) 

( X (3) ,2 (X (2) ,A,E(12,5) ) ) ; 

/* TO PRINT R PRIME */ 

PUT SKIP; 

PUT SKIP EDIT (" REAL PART OF R PRIME= ' , RPR1MER , 

'IMAGINARY PART OF R PRIME- RF'RIMEI ) 

( X (3) ,2 (X (2) ,A,E(12,5) > ) ; 

PUT SKIP; 

PUT SKIP EDIT ( ' ABS (RPR I ME) = ' , RPRABS , ' PHASE (RPR I ME) = 7 , 
RPRPHASE) ( X (3> ,2 ( X (2) , A , E ( 12 , 5) ) ) ; 

PUT SKIP; 

/* TO PTINT T (0) */ 

PUT SKIP; 

PUT SKIP EDIT ( ' REAL PART OF T(0>«',T0R, 

'IMAGINARY PART OF T(0)=',T0I) 

(X (3) ,2 ( X (2) , A , E ( 1 2 , 5 ) ) ) ; 

PUT SKIP ; 

PUT SKIP EDIT < ' REAL PART OF T (-1 ) = ' , TMINUS1R , 

'IMAGINARY PART OF T (— 1 ) = ' , TMINUS1 I ) 

( X ( 3 ) ,2 (X (2) ,A,E(12,5) > ) ; 

END TEMODE; 

/* 

//GO.SYSIN DD * 

. 54 ,.27,. 27, 2. 56, 1.44,1.000,45, 1. 1 
—6 , —5 , —4 , —3 , —2 , — 1 , U , 1 , 2 , o , 4 , 5 
/* 

// 
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